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ABSTRACT 

Phase- shift  keyed  (PSK)  signal  modulation  methods  are 
coming  into  increasing  use  in  modern  communications.  This 
thesis  describes  the  implementation  of  three  methods  of 
computing  the  cyclic  spectrum  to  determine  the  presence  of  PSK 
signals.  The  Strip  Spectral  Correlation  Algorithm  (SSCA)  and 
the  Fast  Fourier  Transform  (FFT)  Accumulation  Method  (FAM) 
both  estimate  the  cyclic  spectral  plane.  The  Sub- FFT 
Accumulation  Method  (SUBFAM)  program  computes  the  Spectral 
Correlation  Function  (SCF)  for  a  set  of  possible  spectral 
frequencies  for  a  single  cycle  frequency.  The  results  of 
algorithm  performance  are  presented  along  with  a  discussion  of 
promising  areas  for  performance  enhancement  and  automation  of 
signal  detection  and  classification. 
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I.   INTRODUCTION 

A .   BACKGROUND 

The  field  of  cyclic  spectral  analysis  is  growing  in 
importance  as  non- stationary  modulation  schemes  proliferate  in 
modern  communications.  Traditionally,  spectral  analysis 
methods  have  assumed  that  signals  of  interest  were  stationary 
or  had  non- time -varying  statistical  parameters.  This 
assumption  is  untrue  for  the  cyclostationary  signals  created 
when  using  pulse  or  carrier  amplitude  modulation,  phase  or 
frequency  carrier  modulation  or  digital  pulse  or  carrier 
modulation.  Gardner  [Ref.  l:pp  419-453]  discusses  many  of 
these  signal  types.  Spread  spectrum  modulation  techniques 
such  as  direct  sequence  PSK  and  frequency- hopped  FSK  also 
exhibit  temporally  varying  spectral  features  [Ref.  l:pp  453- 
457]  .  Reference  2  provides  a  tutorial  on  cyclic  spectral 
analysis . 

Cyclic  spectral  analysis  methods  take  advantage  of 
spectral  variation  over  time  by  correlating  spectral 
estimations  at  discrete  time  intervals  to  locate  signals  that 
may  not  otherwise  be  identified.  A  direct  sequence  signal 
with  a  very  wide  bandwidth  may  be  indistinguishable  from  a 
noisy  background  in  the  frequency  spectral  estimate  of 
Figure  1  [Ref.  3:pp  2-9  -  2-10].    But  this  signal  becomes 


clearly  identifiable  in  the  cyclic  spectral  estimate  of 
Figure  2 . 

One  cyclic  spectral  technique  is  the  frequency  smoothed 
cyclic  periodogram  method  (FSM)  [Ref.  l:pg  464].  This  has 
been  implemented  in  the  Cyclic  Spectral  Analysis  Software 
Package   [Ref.  4]   as  the  program  sxaf.  Hereafter,   this 

software  will  be  referred  to  as  the  SSPI  package  or  the  FSM. 
This  method  was  used  to  generate  Figures  1  and  2 . 
Unfortunately,  the  FSM  is  less  efficient  than  time  smoothing 
methods  for  general  spectral  estimation  [Ref.  5:pg  38]. 

B.   THESIS  GOALS 

The  purpose  of  this  thesis  is  to  implement  two  time 
smoothing  algorithms,  the  FFT  Accumulation  Method  (FAM) 
[Ref.  5:pp  44-47]  and  the  Strip  Spectral  Correlation  Algorithm 
(SSCA)  [Ref.  5:pp  47-48].  The  FAM  is  less  computationally 
complex  than  the  FSM  and  the  SSCA  less  so  than  the  FAM.  Both 
methods  have  been  written  to  be  compatible  with  the  SSPI 
package.  By  integrating  FAM  and  SSCA  into  the  SSPI  package 
future  research  may  more  easily  be  done  in  algorithm 
performance  comparison  and  enhancement. 

Both  the  FAM  and  SSCA  accept  input  data  generated  by  the 
SSPI  utility  program  pam.  They  also  generate  results 

compatible  with  the  SSPI  plotting  utility  plot_sxaf .  Each 
program  may  optionally  generate  output  in  a  format  compatible 
with  the  MATLAB  [Ref.  6]  package  which  is  in  common  use  at  the 


Naval  Postgraduate  School.  Input  and  output  data  are  in  ASCII 
file  formats . 

Throughout  this  thesis  the  same  BPSK  signal  data  is  used 
to  illustrate  the  implementations  of  FAM  and  SSCA.  A 
representative  signal  sample  is  plotted  in  Figure  3.  No  noise 
was  added  to  this  signal.  By  estimating  the  frequency 
spectrum  at  different  points  in  time,  it  is  evident  in  Figure 
4  that  the  spectral  features  do  indeed  vary  temporally. 

Chapter  II  describes  the  implementation  of  the  FAM. 
Chapter  III  describes  the  implementation  of  the  SSCA.  Chapter 
IV  describes  a  utility  to  estimate  cycle  frequencies  for  a 
given  baseband  frequency.  The  Sub-FFT  Accumulation  Method 
(SUBFAM)  program  [Ref.  7]  gives  a  quick  result  for  a  specific 
subset  of  the  spectral  plane.  Chapter  V  discusses  specific 
algorithm  performance  comparisons.  Finally,  Chapter  VI 
presents  conclusions  and  recommendations  for  future  areas  of 
research. 
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Figure  1   Power  Spectral  Density  of  BPSK  Signal  in  Noise 

[Ref.  4:pg  18] 


Figure  2   Cyclic  Spectral  Estimate  of  Direct  Sequence  Signal 
in  Noise  [Ref.  4:pg  18] 
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Figure  3   Typical  BPSK  Signal  Data 


Figure  4   Time  Sequence  Spectrum  of  Typical  BPSK  Signal 


II.   FFT  ACCUMULATION  METHOD 

A .   THEORY 

The  FFT  Accumulation  Method  (FAM)  [Ref.  5:pp  44-47]  is  a 
Fourier  transform  of  correlation  products  between  spectral 
components  smoothed  over  time.  Periodicities  in  the  spectral 
components  then  become  detectable.  The  cyclic  spectral  plane 
as  shown  in  Figure  5  ranges  in  frequency  from  -fs/2  to  +fs/2, 
and  in  cycle  frequency  from  -fs  to  +fs  where  fs  is  the  sampling 
frequency.  For  each  unique  ( f  j ,  afq )  cell  in  Figure  5,  the  FAM 
cross-spectral  estimate  is  defined  to  be  [Ref.  5:pg  44]: 


Sx$Tq     (nL.fi)Lt*Y,XT{xL.fk)Yi{zL,fk)ge{n-x)e      p  (1) 


The   FAM   auto- spectral   estimate   is   defined   to   be 
[Ref.  5:pg  44] : 


P-l  -i2itzq 

siC*  (nL,  fj)  Ac=  £  XT{iL,  fk)  X*(rL,  ft)  gc{n-r)  e       p  (2  ) 


where:  j  is  the  average  frequency  bin  (k+£)/2 

q  is  the  difference  or  bandwidth  frequency  (k--H 
gc(n-r)  represents  the  Hamming  window. 

The  FAM  was  implemented  by  forming  an  array  from  x(kT) 
(0  <  k  <  N-l)  with  rows  which  are  N'  points  long  from  the 
input  sample  data.  The  starting  point  of  each  succeeding  row 
is  offset  from  the  previous  rows  starting  position  by  L 
samples.  A  Hamming  window  is  applied  across  each  row  which  is 
then  Fast  Fourier  transformed  and  downconverted  to  baseband. 
The  result  at  this  point  is  a  two-dimensional  array  with 
columns  representing  constant  frequencies.  Each  column  was 
point-wise  multiplied  in  turn  with  the  conjugate  of  the 
columns  resulting  from  processing  y(kT)  (0  <,  k  <;  N-l)  .  Each 
resultant  product  vector  was  Fast  Fourier  transformed  and  the 
low  frequency  half  placed  into  the  final  cyclic  spectral  plane 
at  the  appropriate  location  in  the  cyclic  spectral  plane. 
Figure  6  shows  a  block  diagram  for  the  FAM  cross- spectral 
estimate . 

The  following  subsections  in  this  chapter  discuss  in 
detail  what  has  been  described  in  the  previous  paragraph.  The 
cycle  frequency  resolution  of  FAM  is  Aa  =  fs/PL  [Ref .  5:pg  44] 
where  fs  is  the  original  data  sampling  rate.  The  frequency 
resolution  of  FAM  is  Af  =  ZjN'  [Ref.  5:pg  44]  .  The  time- 
frequency  resolution  product  is  AtAf  =  PL/N'     [Ref.  5:pg  45]. 


The  command  line  format  for  calling  the  FAM  program  is 
provided  in  Appendix  A.  The  FAM  source  code  listing  is 
provided  in  Appendix  B. 

B .   IMPLEMENTATION 

1.   Input  Channelization 

The  input  sample  data  is  formed  into  a  two-dimensional 
array.  The  array  row  length  is  equal  to  the  number  of  input 
channels  N'  .  For  a  given  number  of  input  sample  points  N;  a 
row  size  of  N'  ,  and  a  chosen  offset  L,  there  are  P  =  (N- 
N'  )  /L  =  N/L  rows  formed.  The  choice  of  N'  must  take  into 
consideration  that  ideally  the  time- frequency  resolution 
product  must  be  much  greater  than  one  [Ref.  5:pg  40].  N' 
should  also  be  a  power  of  2  to  avoid  truncation  or  zero- 
padding  in  the  FFT  routines.  L  should  be  chosen  to  be  less 
than  or  equal  to  N'/4. 

The  completely  filled  array  is  P  rows  by  N'  columns. 
Figure  7  shows  how  a  small  array  is  filled  from  a  discretely 
sampled  signal  x(kT)  when  N' =16 ,  P=8  and  L=4 .  The  number 
inside  each  cell  represents  the  value  of  k  used  to  index  on 
x(kT)  to  fill  that  location  in  the  array. 

Figure  8  shows  the  magnitude  of  the  original  data  for 
the  example  BPSK  signal  organized  into  the  P  by  N'  array  where 
N=512,  N'=32,  L=4  and  P=128.  The  input  data  is  assumed  to  be 
complex  with  a  real  and  imaginary  component  to  each  sample 


point.   Figure  9  shows  a  single  row  of  the  array.   The  phase 
changes  of  the  BPSK  are  evident. 

2 .  Windowing 

A  Hamming  window  [Ref  7:pg  467]  is  applied  to  each  row 
of  the  array.   The  equation  for  the  Hamming  window  is: 

w(r)  =0.54-0.  46cos  (-2?£.)  ,   0  <  r  <  N'-l  (3) 

N'-l 

A  32 -point  Hamming  window  is  plotted  in  Figure  10.  It 
is  applied  to  both  the  real  and  imaginary  parts  of  the  complex 
example  array.  The  magnitude  of  the  resultant  array  is  shown 
in  Figure  11. 

3.  First  FFT 

Each  row  of  the  windowed  data  array  is  Fast  Fourier 
transformed  to  reveal  the  first  spectral  components.  The 
resultant  array  is  still  indexed  P  rows  by  N'  columns  but  now 
the  column  index  relates  to  a  specific  bin  of  spectral 
frequencies.  Figure  12  illustrates  this  relationship. 
Figure  13  shows  the  results  of  FFTing  the  BPSK  example. 

4.  Downconversion 

Each  row  of  spectral  components  is  downconverted  to 
baseband  through  multiplication  with  the  complex  exponential, 

-i2nkmL 

e       N' 


where:   m  is  the  row  index,   0  s  m  <  P-l 

k  is  the  column  index,  0  <,    k  s  N'  - 1 

The  magnitude  of  the  exponential  is  unity  over  the 
array  but  the  phase  shows  considerable  variation.  Figure  14 
shows  the  phase  of  the  exponential  over  the  P  by  N'  array. 
Figure  15  shows  the  phase  for  one  representative  row.  The 
magnitude  of  the  array  remains  unchanged  from  Figure  13 . 

5.  Multiplication 

For  each  cell  in  the  cyclic  spectral  plane,  one  column 
of  the  array  is  multiplied  with  the  conjugate  of  another 
column.  The  separation  of  the  columns  is  determined  by  the 
desired  frequency  separation  or  cycle  frequency  (a=fk-f,) .  The 
mid-point  between  the  columns  is  the  frequency  bin  of  interest 
(f  =  (fk+f i)  /2)  .  Figure  16  shows  two  columns  to  be  conjugate 
multiplied  for  use  in  filling  a  specific  f,a  cell.  Figure  17 
illustrates  the  Z,ot  bin  values  in  the  cyclic  spectral  plane 
for  N'  =4 .  Figure  18  shows  the  corresponding  two  column 
indices  used  to  form  the  product  vector  which  is  used  to 
produce  the  cells  of  Figure  17. 

6 .  Second  FFT 

The  product  from  the  previous  multiplication  is  FFT'd 
to  yield  a  P-point  result.  Only  the  middle  of  the  resulting 
spectrum  is  retained  and  stored  into  the  designated  f ,  a  cell. 
The  upper  and  lower  ends  are  undesireable  because  of  increased 
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estimate  variation  at  the  channel  pair  ends  [Ref.  5:pg  45]. 
Figure  19  illustrates  which  part  of  the  estimate  is  retained 
and  placed  into  the  cell.   Figure  20  shows  the  final  result  of 
the  FAM  auto- spectral  process  on  the  example  BPSK  signal. 
7 .   Data  Reduction 

The  FAM  program  typically  generates  large  output  data 
files.  For  convenience,  an  option  may  be  chosen  to  reduce  the 
amount  of  output.  By  comparison  sorting  for  the  largest  a 
value  in  an  f ,  oi  cell,  the  number  of  at  slices  output  is  reduced 
from  ff'P/2   to  N'  .  Overall   FAM  program  execution  time 

increases  accordingly  to  accomplish  the  sorts. 

Figure  21  illustrates  the  output  full  cyclic  spectral 
plane  storage  array  before  sorting.  Figure  22  shows  the 
output  array  after  data  reduction  has  been  completed.  Figures 
23  and  24  plot  the  resulting  spectral  half-planes  without  and 
with  data  reduction  respectively. 

Since  all  cells  have  an  equal  number  of  a  values,  all 
sorts  are  of  equal  complexity.  Further  work  on  data  reduction 
would  require  selection  of  a  largest  a  value  from  widely 
varying  numbers  of  oi  values  depending  on  the  choice  of  N,  N' 
and  L. 

C .   PERFORMANCE 

An  established  method  of  evaluating  the  complexity  of  an 
algorithm  is  to  determine  the  number  of  floating  point 
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operations  that  must  be  performed.  For  this  estimate  it  is 
assumed  that  an  auto -spectral  estimate  is  being  computed  and 
the  data  is  complex  in  every  step.  It  is  also  assumed  that 
the  Hamming  window  coefficients  and  the  downconversion 
coefficients  have  been  previously  computed  and  stored  for 
later  use.  Each  N-point  FFT  requires  (N/2)*log2N  complex 
floating  point  multiplies  [Ref.  8:pg  506]  or  4* (N/2 ) *log2N 
real  floating  point  multiplies.  The  cost  of  any  output  data 
reduction  is  not  considered  here. 

Applying  the  window 2*P*i\T 

First  FFT 2*P*I\T  *log2N' 

Downconversion 4*P*JV' 

Multiplication 4*N'  *N'  *P 

Second  FFT 2*N'  *Z\T*P*log2P 

Total:        (6+4*i\T  )  *P*N'    +    (2*P*2\T  )  *  (log2W     +   N'  *log2P)  (4) 

Since    P«N/L  (5) 

Total:        (6+4*F)*NF/L    +     (2*NF/L)  (log2F     +    N'  log2N/L)  (6) 
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Figure  5   Generic  ?AM  Cyclic  Spectral  Plane 
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Figure  8   Example  BPSK  Signal  Data  Array,  AT -32,  P=123,  L=4 
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Figure  11   Example  3PSK  Signal  Data  After  Windowing 
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Figure  21   FAM  Storage  Before  Data  Reduction 
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Figure  22   FAM  Storage  After  Data  Reduction 
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Figure  24   FAM  Result  With  Data  Reduction, 
N=512,  N'  =32  and  L=4 . 
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III.   STRIP  SPECTRAL  CORRELATION  ALGORITHM 

A.   THEORY 

The  Strip  Spectral  Correlation  Algorithm  (SSCA) 
[Ref.  5:pp  47-48]  is  a  Fourier  transform  of  correlation 
products  between  spectral  and  temporal  components  smoothed 
over  time.  Periodicities  in  the  spectral  components  then 
become  detectable.  The  cyclic  spectral  plane  as  shown  in 
Figure  25  ranges  in  frequency  from  -fs/2  to  +fs/2,  and  in  cycle 
frequency  from  -fs  to  +£,  where  f,  is  the  sampling  frequency. 
For  each  unique  (fj,aq)  strip  in  Figure  25,  the  SSCA  cross- 
spectral  estimate  is  defined  to  be  [Ref.  5:pg  47] : 


The   SSCA   auto-spectral   estimate   is   defined   to   be 
[Ref.  5:pg  47]  : 


f  a        n~1  -i2nrq 
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where:  k  is  a  multiple  of  the  frequency  resolution, 
-N'  /2  s  k  <  (JV'/2)  -1 
q  is  a  multiple  of  the  cycle  frequency  resolution, 

-N'     <  q  s  (N'  -1) 
g(n-r)  represents  the  Hamming  window. 

The  SSCA  was  implemented  by  forming  an  array  from  x(kT) 
(0  <  k  <  N-l)  with  rows  which  are  N'  points  long  from  the 
input  sample  data.  The  starting  point  of  each  succeeding  row 
is  offset  from  the  previous  rows  starting  position  by  L 
samples.  A  Hamming  window  is  applied  across  each  row  which  is 
then  Fast  Fourier  transformed  and  downconverted  to  baseband. 
The  result  at  this  point  is  a  two-dimensional  array  with 
columns  representing  constant  frequencies.  Each  row  is 
transposed  and  replicated  L  times  for  a  total  of  PL  columns. 
Each  column  is  then  multiplied  by  a  sample  y(kT)  (0  <  k  <  (PL- 
1)  )  .  Each  resultant  row  of  the  array  is  Fast  Fourier 
transformed  and  placed  into  a  strip  of  the  final  cyclic 
spectral  plane  at  the  appropriate  location.  Figure  26  shows 
a  block  diagram  for  the  SSCA  cross -spectral  estimate. 

The  following  subsections  in  this  chapter  discuss  in 
detail  what  has  been  described  in  the  previous  paragraph.  The 
cycle  frequency  resolution  of  SSCA  is  Aa  =  fs/N  [Ref .  5:pg  48] 
where  fs  is  the  original  data  sampling  rate  and  N  is  the 
number  of  input  samples  used.   The  frequency  resolution  of 
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SSCA  is  equal  to  the  sampling  rate  divided  by  the  number  of 
channels  N'  ,  Af  =  tjN'  [Ret.  5:pp  42,48].  The  time- frequency- 
resolution  product  is  AtAf  =  N/N'     [Ret.    5:pg  48] . 

The  command  line  format  for  calling  the  SSCA  program  is 
provided  in  Appendix  C.  The  SSCA  source  code  listing  is 
provided  in  Appendix  D. 

B.   IMPLEMENTATION 

1.   Input  Channelization 

The  input  sample  data  is  formed  into  a  two-dimensional 
array.  The  array  row  length  is  equal  to  the  number  of  input 
channels  N'  .  For  a  given  number  of  input  sample  points  N,  a 
row  size  of  N'  ,  and  a  chosen  offset  L,  there  are  P  =  (N- 
N'  )  /L  "■  N/L  rows  formed.  The  choice  of  N'  must  take  into 
consideration  that  ideally  the  time- frequency  resolution 
product  must  be  much  greater  than  one  [Ref.  5:pg  40].  N' 
should  also  be  a  power  of  2  to  avoid  truncation  or  zero- 
padding  in  the  FFT  routines.  L  should  be  chosen  to  be  less 
than  or  equal  to  N'  /4 . 

The  completely  filled  array  is  P  rows  by  N'  columns. 
Figure  27  shows  how  a  small  array  is  filled  from  a  discretely 
sampled  signal  x(kT)  when  N'  =16 ,  P=8  and  L=4 .  The  number 
inside  each  cell  represents  the  value  of  k  used  to  index  on 
x(kT)  to  fill  that  location  in  the  array. 

Figure  2  8  shows  the  magnitude  of  the  original  data  for 
the  example  BPSK  signal  organized  into  the  P  by  N'   array  where 
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N=512,  N'  =  32,  L=4  and  P=128.  The  input  data  is  assumed  to  be 
complex  with  a  real  and  imaginary  component  to  each  sample 
point.  Figure  29  shows  a  single  row  of  the  array.  The  phase 
changes  of  the  BPSK  are  evident. 

2 .  Windowing 

A  Hamming  window  [Ref  7:pg  467]  is  applied  to  each  row 
of  the  array.   The  equation  for  the  Hamming  window  is: 

w(r)  =0  .  54-0  .46 cos  (  2nr  )  ,   0  <L  r  <  N'-l  (9) 

N'-l 

A  32 -point  Hamming  window  is  plotted  in  Figure  30.  It 
is  applied  to  both  the  real  and  imaginary  parts  of  the  complex 
example  array.  The  magnitude  of  the  resultant  array  is  shown 
in  Figure  31. 

3.  First  FFT 

Each  row  of  the  windowed  data  array  is  Fast  Fourier 
transformed  to  reveal  the  first  spectral  components.  The 
resultant  array  is  still  indexed  P  rows  by  N'  columns  but  now 
the  column  index  relates  to  a  specific  bin  of  spectral 
frequencies.  Figure  32  illustrates  this  relationship. 
Figure  33  shows  the  results  of  FFTing  the  BPSK  example. 

4.  Downconversion 

Each  row  of  spectral  components  is  downconverted  to 
baseband  through  multiplication  with  the  complex  exponential, 
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where:   m  is  the  row  index,   0  s  m  s  P-l 

k  is  the  column  index,  0  <.   k  s  N'  - 1 

The  magnitude  of  the  exponential  is  unity  over  the 
array  but  the  phase  shows  considerable  variation.  Figure  34 
shows  the  phase  of  the  exponential  over  the  P  by  N'  array. 
Figure  35  shows  the  phase  for  one  representative  row.  The 
magnitude  of  the  array  remains  unchanged  from  Figure  33. 

5.  Replication 

Each  row  is  copied  into  one  column  of  an  empty  N'  by 
PL  array.  It  is  then  replicated  in  the  L-l  adjacent  columns. 
Figure  36  illustrates  this  process. 

6.  Multiplication 

Each  column  of  the  array  is  pointwise  multiplied  with 
the  conjugate  of  a  sample  value  y(kT) .  There  are  PL=N  columns 
and  PL  samples  from  y(kT).  Figure  37  illustrates  the 
conjugate  multiplication  process. 

7 .  Second  FFT 

Each  row  from  the  previous  multiplication  is  Fast 
Fourier  transformed  to  yield  a  PL-point  result.  Each 
resulting  vector  is  stored  into  a  strip  of  the  cyclic  spectral 
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plane.   Figure  3  8  shows  the  final  result  of  the  SSCA  auto- 
spectral  process  on  the  example  BPSK  signal. 
8.   Data  Reduction 

The  SSCA  Program  typically  generates  large  output  data 
files.  For  convenience,  an  option  may  be  chosen  to  reduce  the 
amount  of  output.  By  comparison  sorting  for  the  largest  a 
value  in  an  f,oi  cell,  the  number  of  a  slices  is  reduced  from 
N  to  N/L.  Overall  SSCA  execution  time  increases  accordingly 
to  accomplish  the  searches. 

Figure  39  illustrates  the  output  full  cyclic  spectral 
plane  storage  array  before  sorting.  Figure  40  shows  the 
output  array  after  data  reduction  has  been  completed.  Figures 
41  and  42  plot  the  resulting  spectral  half-planes  without  and 
with  data  reduction  respectively. 

Since  all  cells  have  an  equal  number  of  a  values,  all 
sorts  are  of  equal  complexity.  Further  work  on  data  reduction 
would  require  selection  of  a  largest  oi  value  from  widely 
varying  numbers  of  a  values  depending  on  the  choice  of  N,  N' 
and  L. 

C .   PERFORMANCE 

An  established  method  of  evaluating  the  complexity  of  an 
algorithm  is  to  determine  the  number  of  floating  point 
operations  that  must  be  performed.  For  this  estimate  it  is 
assumed  that  an  auto- spectral  estimate  is  being  computed  and 
the  data  is  complex  in  every  step.   It  is  also  assumed  that 
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the  Hamming  window  coefficients  and  the  downconversion 
coefficients  have  been  previously  computed  and  stored  for 
later  use.  Each  N-point  FFT  requires  (N/2)*log2N  complex 
floating  point  multiplies  [Ref.  8:pg  506]  or  2*N*log2N  real 
floating  point  multiplies.  The  cost  of  any  output  data 
reduction  is  not  considered  here. 

Applying  the  window 2*P*N' 

First  FFT 2*P*I\T  *log2^' 

Downconversion 4*P*W 

Multiplication 4*N*i\T 

Second  FFT 2*N*I\T  *log:N 

Total:       2*N'*(6*P    +    4*N    +     (2*P    +    2*N)*log2N)  (10) 

since      P«N/L  (11) 

Total:       2*i\T  *  (  (6*N/L    +    4*N)     +     (2*N/L    +    2*N)*log2N)  (12) 
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Figure   27   Layout   of  A  Sample   Input  Array  Showing   Input 
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Figure  28   Example  BPSK  Signal  Data  Array,  AT  =32,  P=12  3,  L=4 
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Figure  31   Example  BPSK  Signal  Data  After  Windowing 
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Figure  32   Generic  Array  Row  After  First  FFT 
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Figure  34  Phase  of  the  Downconversion  Exponential 
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Figure  3  8   Final  SSCA  Results  of  the  BPSK  Example, 
N=512,  N'  =32  and  L=4. 
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Figure  39   SSCA  Storage  Array  Before  Data  Reduction 
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Figure  40   SSCA  Storage  Array  After  Data  Reduction 
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Figure  41   BPSK  Example  Results  Without  Data  Reduction, 
N=512,  N'  =32  and  L=4 . 
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Figure  42   BPSK  Example  Results  With  Data  Reduction, 
N=512,  N'  =32  and  L-4 . 
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IV.   SUB-FFT  ACCUMULATION  METHOD 

A.   THEORY 

In  the  course  of  this  thesis  work  a  need  arose  for  a 
flexible  and  customized  Spectral  Correlation  Function  (SCF) 
[Ref.  7]  to  determine  the  presence  of  phase-shift  keyed 
signals  in  experimental  data.  The  SCF  [Ref.  7]  is 
represented: 

N-l  .l2n  ll/  ,     2nr  , 

sZy(£)=^x{r)y{i)e        r*    OSFxOSFve^     N  (13) 

1  =  0 

where : 

N  is  the  number  of  data  samples  used 

f,  is  the  intermediate  modulation  frequency 

fs  is  the  sampling  frequency 

OSF  is  a  one-sided  bandpass  filter 

The  program  written  to  implement  this  algorithm  is  named 
the  Sub-FFT  Accumulation  Method  (SUBFAM)  program  to 
distinguish  it  from  the  more  general  purpose  FAM  program.  The 
command  line  format  for  calling  the  SUBFAM  program  is  provided 
in  Appendix  E.  The  SUBFAM  source  code  is  provided  in 
Appendix  F. 
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B .  IMPLEMENTATION 

Figure  43  illustrates  the  operations  performed  in  the 
SUBFAM  program.  A  set  of  N  points  is  obtained  from  within  a 
file  of  signal  data  values.  Each  set  of  N  points  is  processed 
through  a  one-sided  bandpass  filter  to  remove  either  all 
positive  or  negative  frequencies  as  desired.  Downconversion 
of  each  file  to  baseband  from  the  original  sampling  and 
transmission  frequency  bands  follows.  The  results  are  then 
correlated  through  a  complex  conjugate  multiplication.  A 
final  N-point  Fast  Fourier  Transform  yields  the  spectral 
correlation  result. 

C .  PERFORMANCE 

The  SUBFAM  program  performed  as  required.  Figure  44 
illustrates  the  results  of  processing  test  data  containing 
phase- shift -keyed  signal  samples.  The  peaks  at  the  chip 
frequencies  are  apparent.  Figure  45  shows  the  results  of 
correlating  test  data  with  data  containing  only  white  noise. 
The  lack  of  correlation  between  spectral  features  yields 
results  which  are  four  orders  of  magnitude  less  than  in  the 
successful  signal  detection  of  Figure  44. 
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Figure    43      SUBFAM  Program  Block  Diagram 
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Figure   44      SUBFAM  Program  Result  With  Signal,    N=409S 
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Figure   45      SUBFAM  Program  Result  With  Noise,    N=4096 
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V.   ALGORITHM  PERFORMANCE 

The  complexity  of  the  FAM  and  SSCA  algorithms  is  estimated 
in  Chapter  II  Section  C  and  Chapter  III  Section  C 
respectively.  Both  require  on  order  of  Nlog:N  floating  point 
multiplies  where  N  is  the  original  number  of  sample  points 
used.  Brown  and  Loomis  [Ref .  9]  explore  the  complexity  of  the 
FAM,  SSCA  and  FSM  algorithms  in  some  detail.  Their  results 
are  consistent  with  an  order  Nlog:N  multiplies  for  FAM  and 
SSCA. 

The  cyclic  spectral  plane  is  symmetric  about  the  f=0  and 
the  oi=Q  axes  for  real  signals.  For  the  most  efficient 
performance  it  is  only  necessary  to  compute  one  quadrant  of 
the  cyclic  spectral  plane.  To  compare  the  performance  of  FAM 
and  SSCA  with  FSM,  complexity  figures  derived  from  Reference 
9  for  quadrant  complexity  will  be  used.  The  three  algorithms 
require  the  following  numbers  of  floating  point  multiplies: 

FAM  Complexity  [Ref.  9:pg  18] 

2Ntf'log2(  — )  +  8Nlog.  (N')    +  ANN'   =  20N       (14) 
N' 
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SSCA   Complexity     [Ref.    9 : pg    20] 

N^'log2(N)    +  2N7V7  +  SNlog^iV7)    +  12N  (15) 

FSM  Complexity    [Ref.    9:pg    16] 

N2   +  2Nlog2  (N)  (16) 

The  complexity  of  four  example  cases  with  varying  N  and  N' 
are  estimated  below  using  equation  (14),  (15)  and  (16) . 

FAM SSCA FSM 


1024 

64 

5.3*105 

3  .6*105 

106 

2048 

128 

2.0*106 

1.4*106 

4.2*106 

32, 768 

512 

1.5*108 

1.1*108 

109 

1, 048, 

576 

8192 

8.0*1010 

7.0*1010 

1012 

It  is  evident  that  FAM  and  SSCA  require  substantially 
fewer  multiplies  than  FSM.  Particularly  with  larger  values  of 
N,  FAM  and  SSCA  require  Nlog:N  while  FSM  requires  N2  floating 
point  multiplies. 
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VI.   CONCLUSIONS 

A.  SUMMARY 

The  purpose  of  this  thesis  was  to  implement  the  FAM  and 
SSCA  for  use  in  cyclic  spectral  analysis  research.  Both 
algorithms  were  successfully  implemented  in  a  manner 
compatible  with  the  SSPI  analysis  package  to  facilitate  future 
work.  By  allowing  the  user  a  choice  of  two  output  file 
formats,  results  may  easily  be  used  in  either  MATLAB  [Ref.  6] 
or  SSPI  [Ref.  4]  for  further  signal  processing  applications. 

It  has  been  shown  that  the  FAM  and  SSCA  algorithms 
generate  results  consistent  with  FSM  but  require  substantially 
fewer  floating  point  multiplies  than  FSM.  This  is  especially 
true  as  the  number  of  sample  data  points  used  increases. 

B .  RECOMMENDATIONS 

In  the  course  of  this  thesis  it  has  become  apparent  that 
the  FAM  and  SSCA  algorithms  consist  of  inherently  parallel 
operations.  Roberts  and  Loomis  explore  these  possibilities  in 
Reference  10.  The  high-level  sequential  nature  of  these 
algorithms  also  lend  them  to  pipelining  architectures.  By 
combining  a  parallel  or  multi -processor  approach  with 
pipelining  a  large  improvement  in  performance  should  be 
obtainable.   This  would  be  espectially  true  in  applications 
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requiring  large  sets  of  signal  data  samples  or  repeated, 
time -sequenced  cyclic  spectral  plane  estimations. 

Two  areas  which  could  benefit  from  further  work  are 
automatic  signal  identification  and  signal  parameterization 
algorithms.  The  results  from  cyclic  spectral  analysis 
algorithms  such  as  FAM  and  SSCA  may  be  input  into  Linear 
Associator  or  Mapping  neurocomputer  networks  as  described  in 
Reference  11.  Intensity  transformations  using  digital  image 
processing  techniques  from  Reference  12  also  appear  to  have 
promise.  Both  classes  of  techniques  could  have  utility  in 
automatic  signal  identification  and  parameterization.  They 
also  have  the  added  benefit  that  they  lend  themselves  to 
parallel  and  pipelined  computing  architectures. 
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APPENDIX  A.   FAM  PROGRAM  USE 


Correct  commands  are: 

for  the  cross- fam 

fam  inputfilel  inputfile2  outputfile  n  nprime  1  Oflag 
fam  inputfilel  inputfile2  outputfile  n  nprime  1  Oflag  red 

or  for  the  auto- fam 

fam  inputfilel  outputfile  n  nprime  1  Oflag 
fam  inputfilel  outputfile  n  nprime  1  Oflag  red 

where:  inputfilel  is  the  file  containing  the  x  signal  samples 
inputfile2  is  the  file  containing  the  y  signal  samples 
outputfile  is  where  the  spectrum  values  will  be  placed 
n  is  the  number  of  samples  to  use  from  the  inputfiles 

and  it  must  be  a  power  of  2 
nprime  is  the  group  size  of  the  input  datasets  and  it 

must  be  a  power  of  2 
1  is  the  starting  offset  of  subsequent  datasets  and 

it  must  be  a  power  of  2 
Oflag  indicates  the  output  filetype,  ascii  or  binary 
-asc  for  an  ascii,  MATLAB  compatible  file, 

full  cyclic  spectral  plane 
-plo  for  a  plot_sxaf,    SSPI  compatible  file, 
half  cyclic  spectral  plane 
red  reduces  the  amount  of  data  output 

note:   the  first  two  entries  in  every  input  file  are  expected 
to  be  the  datatype  and  n.   datatype  =  1  for  real 
values,  2  for  complex  values.   n  is  the  number  of 
samples  contained  in  the  file,  one  sample  per  line. 

input  file  format: 

datatype    n 
sample  #1 
sample  #2 


sample  #n 
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output  file  ASCII  MATLAB  format 

nprime      n 
output (1)  (1) 
output (1)  (2) 


output (nprime) (n) 


output  file  SSPI  plot_sxaf   format: 

datatype  (1  for  real,  2  for  complex) 
number_of_alphasalpha_minalpha_max 
number_of_f reqsf req_minf req_max 

value  of  alpha  #1 

number  of  freqs  in  #lf req_minf req_max 

spectrum  at  freq_min 


spectrum  at  freq_max 


value  of  alpha  #alpha_max 

number  of  freqs  in  #alpha_maxf req_minf req_max 

spectrum  at  freq_min 


spectrum  at  freq_max 
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APPENDIX  B.   FAM  PROGRAM  LISTING 
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^include  <stdlib.h> 
•tmclude  <stdio.h> 
^include  <math.h> 

# include  "/home3/carter/thesis/SSPI/pam/f ft . c" 
# include  " /home 3 /carter /thesis /SSP I /pam/ radix, c" 
main (argc, argv) 
int  argc; 
char  *argv [ ] ; 
{ 

COMPLEX  *x, **s3, **y3,  **dwnconv,  **s4; 
COMPLEX  *tempf ft, *yl, *sl, **s2, **y2, * window ; 
/*   x  passes  data  to/from  fft  routine 
s3  holds  results  of  FFTmg  s2 
y3  holds  results  of  FFTmg  y2 
dwnconv  holds  downconversion  coefficients 
s4  holds  correlation  multiply  results 
tempfft  passes  data  to/from  fft  routine 
yl  holds  y  samples  from  mputfile  2 
si  holds  x  samples  from  inputfile  1 
s2  holds  channelized  x  samples 
y2  holds  channelized  y  samples 
window  holds  Hamming  window  values 

*/ 
int  i, j,k, type, n, file_n, nprime,p, 1, direction,  norm; 
int  a,b,  c,  cross, num_f,max_num_f, data_type,max _num_alf ; 
int  tempi, temp], half p, reduce, curr_alf ; 
int  *outint; 

float  numreal,  numimag, convfac,mainfac, num_alf , f_min, f_max; 
float  alpha_min,  alpha_max,  f_mm_all,  f_max_all; 
float  tempreal, tempimag, bigmag, tempmag; 
float  twopi=6. 28318530718; 
float  *outreal; 
double  z,y; 

char  *mfilel,  *mfile2,  *outfile; 
FILE  *ifpl, *ifp2, *ofp; 

/*  check,  for  the  correct  number  of  input  arguments 
the  correct  commands  are: 

fam  mputfilel  mputfile2  outputfile  n  nprime  1  Oflag  reduce 

or 
fam  mputfilel  inputfile2  outputfile  n  nprime  1  Oflag 

or 
fam  mputfilel  outputfile  n  nprime  1  Oflag  reduce 

or 
fam  mputfilel  outputfile  n  nprime  1  Oflag 

where:   mputfilel  is  where  the  x  signal  samples  are  expected  to  be 
inputfile2  is  where  the  y  signal  samples  are  expected  to  be 
outputfile  is  where  the  spectrum  values  will  be  put 
n  is  the  number  of  samples  to  use  from  the  mputfiles 
nprime  is  the  group  size  of  input  datasets 
1  is  the  starting  offset  of  subsequent  datasets 
Oflag  indicates  the  output  filetype,  ascii  or  binary 

-asc   for  an  ascii,  matlab  compatible  file,  full  plane 
-plo   for  a  plot_sxaf  compatible  file,  half  plane 
reduce  is  an  option  to  reduce  the  amount  of  output 
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iatatype  and  n,  are  expected  to  ce  at  the 

datatype  =  1  for  real,  2  for  complex 
n  is  the  number  of  samples  foilowincr 


>p  of  the  incut 


.le 


written  by  LCDR  Nancy  J.  Carter,  USN    NPGS  Monterey  CA   01  October  1992 
to  implement  the  Frequency  Accumulation  Method  of  cyclic  spectral  analysis 
as  described  in  the  Brown/Gardner/Loomis  IEEE  paper 


20  October  1992  ....  changed  to  ascn  MAT  LAB  output  format 

22  October  1992 ....  added  plot_sxaf  output  format  compatible  with  SSPI 

23  November  1 992 ....  added  option  to  reduce  amount  of  output 

*/ 
CHECK  INPUT  VALUES     */ 


Oflagxn" ) 
); 


9))  { 
'fam. ...  fatal  errorXn"); 

incorrect  number  of  calling  arguments \n" ) ; 

correct  formats  are:\n"); 

...   fam  mputfilel  inputfile2  outputfile  n  nprime  1  Oflag  reduce\n") 
fam  mputfilel  mputfile2  outputfile  n  nprime 
fam  inputfile  outputfile  n  nprime  1  Oflag  reduceXn' 

...   fam  inputfile  outputfile  n  nprime  1  OflagXn"); 
where \n" ) ; 

...      mputfilel  contains  the  x  samples \n" ) ; 

...      mputfile2  contains  the  y  samplesXn"); 

...      outputfile  will  contain  the  resultsXn") ; 

...       n  is  the  number  of  input  samples  to  useXn"); 

...      nprime  is  the  group  size  of  input  datasetsXn" ) ; 

...      1  is  the  offset  of  subsequent  datasets \n" ) ; 

...      Oflag  is  the  output  file  formatXn"); 

...  Oflag  =  -asc,  an  ascn  file  is  producedXn" ) ; 

...  Oflag  =  -plo,  a  plot_sxaf  file  is  producedXn"); 

...      reduce  is  an  option  for  reduced  amount  of  outputXn"); 


/' 

/*  */ 

if  (  (argc  !=  7)&&(argc  !=  8)<£&(argc 

printf ( 

prmtf  ( 

printf ( 

printf ( 

printf ( 

printf ( 

printf ( 

printf ( 

printf ( 

printf ( 

printf ( 

printf ( 

printf ( 

printf ( 

printf ( 

printf ( 

printf ( 

printf ( 

exit  ()  ; 

} 
if  (argc==7)  {   /' 

cross=0; 

inf ilel=argv [  1  ]  ; 

outfile=argv[2] ; 

n=atoi (argv[3] ) ; 

nprime=atoi (argv [4]  )  ; 

l=atoi (argv[5] )  ; 

reduce=0; 

} 
if  ( (argc==8)&& (argv [7] [0] !=' r' ) ) 

cross=l; 

inf ilel=argv [ 1 ] ; 

mfile2=argv[2]  ; 

outf ile=argv [3] ; 

n=atoi (argv [4] )  ; 

nprime=atoi (argv [ 5 ] ) ; 

l=atoi (argv [6] ) ; 

reduce=0; 
} 
if  ( (argc==8)&& (argv [7] [0]=='r' ) ) 

cross=0; 


this  is  an  auto-fam,  no  output  reduction 


/*  this  is  a  cross-fam,  no  output  reduction 


/*  this  is  an  auto-fam,  with  output  reduction 
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mf  ilel=argv  [1] 
outf ile=argv [2] 
n=atoi (argv[3] ) 
nprime=atoi  (arg- 
l=atoi (argv[5] ) 
reduce=l; 


v[4]); 


if  (argc==9)  {   /*  this  is  a  cross-fam,  with  output  reduction   * 

cross=l; 

infilel=argv [1] ; 

mfile2=argv[2]  ; 

outfile=argv[3] ; 

n=atoi (argv[4] ) ; 

nprime=atoi (argv [ 5 ] ) ; 

l=atoi (argv [6] ) ; 

reduce=l; 

} 
/*  verify  that  n,  nprime  and  1  are  powers  of  2      */ 
i=n%2; 
if  (i!=0) { 

printf ("fam. . . . fatal  error\n" ) ; 

printf  (" calling  argument  n  is  not  a  power  of  2\n"); 

exit  () ; 

} 
i=npnme%2; 
if  (i!=0){ 

printf ("fam. ...  fatal  error\n"); 

printf  (" calling  argument  nprime  is  not  a  power  of  2\n") 

exit  ()  ; 

} 
j-1%2; 
if  (j'.-OM 

printf ("fam. . . . fatal  error\n" ) ; 

printf  (" calling  argument  1  is  not  a  power  of  2\n"); 

exit  ()  ; 

} 
/*  */ 

/*  open  input  file  1,  prepare  to  get  the  x  signal  sample  data 
/*  */ 

lfpl  =  fopen (infilel, "r") ; 
fscanf (ifpl, "%i  %i", &type, &file_n) ; 

/*  verify  that  file_n  is  greater  than  or  egual  to  n  */ 
if  (file_n<n) { 

printf ("fam. ...  fatal  error\n"); 

printf (" inputfilel  does  not  contain  enough  samples\n"); 

fclose (ifpl) ; 

exit () ; 

} 
/*  find  p  -  the  number  of  datasets (nprime  long)   */ 
p= ( (n-nprime) /l) ; 

/*  */ 

/*     */ 

/*     READ  IN  DATA  SAMPLES    */ 
/*     */ 

/*  allocate  space  to  read  in  the  sample  values  */ 
/*     */ 
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sl= (COMPLEX*) calloc(n/  sizeof (COMPLEX) )  ; 
if  (sl==NULL) { 

prmtf  ("fam.  ...  insufficient  space  to  allocate  si\n"); 
exit  ( ) ; 
} 
if  (type— 1)  { 

/*  */ 

/*   read  in  the  real  sample  values  */ 

/*  */ 

for  (i=0;  l  <  n;  i++) { 

f scanf (ifpl, "%e\n", snumreal) ; 
si [i] . r=numreal; 
si  [i]  .1=0.0; 
} 
} 

else  { 

/*  */ 

/*   read  in  the  complex  sample  values  */ 

/*  */ 

for  (i=0;  i  <  n;  i++)  { 

fscanf (ifpl, "%e   %e\n" , inumreal, inumimag) ; 
si [i] .r=numreal; 
si [i] . i=numimag; 
} 
} 

/*     */ 

/*   close  the  input  file  #1    */ 
/*     */ 
fclose (ifpl) ; 
I*  */ 

/*  if  this  is  a  cross-fam  go  get  y  samples   */ 
/*     */ 
if  (cross==l)  { 

/*  V 

/*  open  mputfile2,  prepare  to  get  the  y  signal  sample  data   */ 
/*  */ 

ifp2  =  fopen (inf ile2, "r") ; 
fscanf  (ifp2, "%i  %i", Stype, &file_n) ; 

/*  verify  that  file_n  is  greater  than  or  equal  to  n  */ 
if  (file_n<n) { 

prmtf  ("fam.  .  .  .  fatal  error\n"  )  ; 

print f (" input filel  does  not  contain  enough  samples\n" ) 

fclose  (ifpl) ; 
exit  ()  ; 
} 
/*     */ 

/*     READ  IN  Y  DATA  SAMPLES    */ 
/*     */ 

/*  allocate  space  to  hold  the  sample  values  */ 
/*     */ 

yl= (COMPLEX*) calloc(n, sizeof (COMPLEX) )  ; 
if  (yl==NULL) { 

prmtf  ("fam.  ...  insufficient  space  to  allocate  yl\n"); 
exit  ()  ; 
} 
/*  */ 
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/*   read  in  the  complex  sample  values  */ 

/*  */ 

for  (i=0;  i  <  n;  i++)  { 

fscanf (ifp2, "%e   %e\n", snumreal, inumimag)  ; 
yl[i] . r=numreal; 
y 1 [ i ] . i=numimag ; 
} 
/*     */ 

/*   close  the  input file2    */ 
/*     */ 

fclose(ifp2) ; 
}   /*  end  if  cross=l    */ 
/*     */ 

/*     FORM  DATASETS     */ 
/*     */ 

/*  allocate  space  to  hold  the  p-by-nprime  datasets   */ 
/*     */ 

s2= (COMPLEX**) calloc (p, sizeof (COMPLEX*) ) ; 
if  (s2==NULL) { 

prmtf  ("fam insufficient  space  to  allocate  s2\n"); 

exit  ()  ; 
} 
for  (i=0;  i<p;  i++)  { 

s2 [i]  =  (COMPLEX*) calloc  (nprime, sizeof (COMPLEX) )  ; 
if  (s2[i]==NULL) { 

printf ("fam. ...  insufficient  space  to  allocate  s2\n"); 
exit  (); 
} 
} 
/*  */ 

/*  form  p  datasets  of  nprime  samples  from  the  original  sample  stream   */ 
/*  */ 

for  (i=0;  i<p;  i++)  { 

for  (j=0;  ]<nprime;  j++) { 
k=i*l+j; 
S2[i]  [j]=sl[k]; 
} 
} 

/*     */ 

/*  if  this  is  a  cross-fam  form  y  datasets  */ 
/*     */ 
if  (cross==l)  { 

/*  */ 

/*     */ 

/*  allocate  space  to  hold  the  p-by-nprime  datasets   */ 
/*    */ 

y2= (COMPLEX* *) calloc (p, sizeof (COMPLEX*) ) ; 
if  (y2==NULL) { 

printf ("fam. ...  insufficient  space  to  allocate  y2\n"); 
exit  ()  ; 
} 
for  (i=0;  Kp;  i++)  { 

y2 [i]= (COMPLEX*) calloc (nprime, sizeof (COMPLEX) ) ; 
if  (y2 [i]==NULL) { 

printf ("fam. ...  insufficient  space  to  allocate  y2\n"); 
exit ( ) ; 
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} 
} 
/*  */ 

/*  form  p  datasets  of  nprime  samples  from  the  original  sample  stream   */ 

*/ 
for  (i=0;  i<p;  i++)  { 

for  (j=0;  ]<nprime;  ]++){ 
k=i*l+i; 
y2[i]  [j]-yl[k]; 
} 
} 
}   /*  end  if  cross=l    */ 
/*    V 

/*     APPLY  WINDOW  TO  DATASETS     *  I 
/*     */ 

/*  allocate  space  to  hold  the  window  multiplicand  values  */ 
/*     */ 

wmdow=( COMPLEX*)  calloc  (nprime,  sizeof  (COMPLEX)  )  ; 
if  (wmdow==NULL)  { 

prmtf  ("fam insufficient  space  to  allocate  sl\n"); 

exit  ()  ; 
} 
/*  V 

/*  calculate  the  window  values  for  the  "nprime"  sample  wide  window   */ 
/*  */ 

for  (i=0;  Knprime;  i++)  { 
y= (twopi*i) / (nprime-1) ; 
z=cos (y) ; 

window [l] .r=  0.54  -  (0.46*z); 
window [i].i=  window [i].r ; 
} 

/*  */ 

/*  apply  the  window  to  all  rows  of  s2       */ 
/*  */ 

for  (i=0;  Kp;  i++)  { 

for  (j=0;  ]<nprime;  j++) { 

s2 [i]  [ j]  .r=s2  [i]  [  j]  .r*window[]]  .  r; 
s2 [i] [ j] . i=s2 [i] [ j] . i 'window [j] . i; 
} 

;.   ./ 

/*  if  this  is  a  cross-fam  apply  window  to  y2  */ 
/*     */ 
if  (cross==l)  { 

/*  */ 

for  (i=0;  i<p;  i++) { 

for  (j=0;  ]<nprime;  j++) { 

y2 [i]  []]  .r=y2  [l]  [  j]  .r*window[ j]  .r; 
y2[i] [j] .i=y2[i] [j] . i*window[ j ] .i; 
} 
} 
}   /*  end  if  cross=l    */ 
/*     */ 

/*   FFT  EACH  DATASET  ROW     */ 
/*     */ 
/*   allocate  space  to  hold  rows  of  data  for  passing  to  the  FFT  routine  */ 
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*/ 
x= (COMPLEX*) calloc (nprime, sizeof (COMPLEX) ) ; 
if  (x==NULL) { 

prmtf  ("fam insufficient  space  to  allocate  x\n"); 

exit  ()  ; 

} 
/*     */ 

/*   allocate  space  to  hold  rows  of  results  from  the  FFT  routine  */ 
/*     */ 

s3= (COMPLEX**) calloc (p, sizeof (COMPLEX*)  )  ; 
if  (s3==NULL) { 

prmtf ("fam. ...  insufficient  space  to  allocate  s3\n"); 

exit  ( ) ; 

} 
for  (i=0;  i<p;  i++)  { 

S3 [l]  =  (COMPLEX*) calloc (nprime, sizeof (COMPLEX) )  ; 

if  (s3 [i]==NULL) { 

prmtf  ("fam.  ...  insufficient  space  to  allocate  s3\n"); 

exit  ( ) ; 

} 

} 
/*  */ 

/*  take  an  FFT  of  each  row  of  s2  */ 

/*  */ 

/*    set  values  in  arguments  sent  to  fft   */ 
direct ion=l; 
norm=l; 
for  (i=0;  i<p;  i++)  { 

/*    copy  row  of  s2  into  complex  array      */ 

for  ( j=0;  ]<nprime;  j++) { 
x[j]=s2[i]  []]; 
\ 

/*    go  get  FFT  performed  */ 

fft (x, nprime, direction, norm) ; 

/*  copy  x  values  into  a  row  of  s3  */ 

for  (j=0;  ]<nprime;  j++) { 
s3[i] [j]=x[j]; 
} 
} 

/*     */ 

/*  if  this  is  a  cross-fam  take  an  FFT  of  each  row  of  y2  */ 
/*    */ 
if  (cross==l)  { 

/*  */ 

/*     */ 

/*   allocate  space  to  hold  rows  of  results  from  the  FFT  routine  */ 
/*     */ 

y3= (COMPLEX**) calloc (p, sizeof (COMPLEX*) ) ; 
if  (y3==NULL) { 

prmtf  ("fam.  ...  insufficient  space  to  allocate  y3\n"); 

exit  ()  ; 

} 
for  (i=0;  i<p;  i++) { 

y3 [i]  =  (COMPLEX*) calloc (nprime, sizeof (COMPLEX) )  ; 

if  (y3[i]==NULL) { 

prmtf  ("fam.  ..  .insufficient  space  to  allocate  y3\n"); 
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exit  ()  ; 


for  (i=0;  i<p;  i++)  { 
/*     copy  row  of  y2  into  complex  array      */ 
for  (j=0;  j<nprime;  j++) { 
x[j]-y2[i]  [j]; 
} 
/*     go  get  FFT  performed  */ 

f  ft  (x,  npnme,  direction,  norm)  ; 

/*  copy  x  values  into  a  row  of  y3  */ 

for  (j=0;  ]<nprime;  j++)  { 
y3[i]  t]]=x[]j; 
} 
} 
}   /*  end  if  cross=l    */ 
/*     */ 
/*     */ 

/*     DOWNCONVERT  EACH  ROW     */ 
/*      */ 

/*   allocate  space  to  hold  the  downconversion  multiplicands   */ 
/*     */ 

dwnconv= (COMPLEX**) calloc (p, sizeof (COMPLEX* ) ) ; 
if  (dwnconv==NULL) { 

printf ("fam. ...  insufficient  space  to  allocate  dwnconv\n"); 
exit  ()  ; 
} 
for  (i=0;  i<p;  i++)  { 

dwnconv  [i]  =  (COMPLEX*)  calloc  (npnme,  sizeof  (COMPLEX)  )  ; 
if  (dwnconv [l] ==NULL) { 

printf ("fam. ...  insufficient  space  to  allocate  dwnconv\n"); 
exit  ()  ; 
} 
} 
/*  */ 

/*   downconvert  each  of  the  transform  sequences     */ 
/*  */ 

/*      calculate  the  down  conversion  multipliers     */ 
/*  */ 

mainfac=twopi*l/nprime; 
for  (i=0;  i<p;  i++) { 

for  ( j=0;  ]<nprime;  j++) { 
convfac=i* j*mainfac; 
dwnconv [i]  [ j ]  . r=cos (convfac)  ; 
dwnconv [i] [ j ] . i=  (-1 . 0) *sin (convfac) ; 
} 
) 

/*  multiply  downconversion  factor  against  each  frequency  value   */ 
for  (i=0;  i<p;  i++)  { 
for  (]=0;  ]<nprime;  j++) { 

numreal=s3 [i]  [ j ]  . r*dwnconv [i]  [j].r  -  s3  [i] [ j] . i*dwnconv[i]  [ j] . i; 
s3 [i]  [ j ]  . i=s3 [i]  [ j ]  .r*dwnconv[i]  [j].i  +  s3  [i]  [ j] . i*dwnconv [i]  [ j]  . r; 
s3 [i] [ j ] . r=numreal; 
\ 
\ 
if  (cross==l)  { 
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*/ 

/*  if  this  is  a  cross-fam  downconvert  v3  */ 
/*     */ 

/*  multiply  downconversion  factor  against  each  frequency  value   */ 
for  (i=0;  i<p;  i++)  { 

for  (j=0;  j<npnme;  j++)  { 

numreal=y3  [i]  [j]  . r*dwnconv [i ]  [j]  .r   -   y3[i]  [j]  . i*dwnconv [i]  [j]  .1; 
y3[i] [j] .i=y3[i] [j] . r*dwnconv [i] [j] .1  +  y3[i] [j] . i*dwnconv [i] [j] .r; 
y3 [i] [ j ] . r=numreal; 
} 
} 
}      /*  end  if  cross=l   */ 
/*     */ 

/*     MULTIPLY  COLUMNS     */ 
/*     */ 

/*   allocate  space  to  hold  the  correlation  multiply  (column  multiply)  results  */ 
/*     */ 
a=nprime*p/2; 

s4= (COMPLEX**) calloc (a, sizeof (COMPLEX*) ) ; 
if  (s4==NULL) { 

printf (" fam.  ...  insufficient  space  to  allocate  s4\n"); 
exit  ()  ; 
} 
for  (i=0;  Ka;  i++)  { 

s4 [i] = (COMPLEX*) calloc (nprime, sizeof (COMPLEX) ) ; 
if  (s4 [i]==NULL) { 

printf ("fam. ...  insufficient  space  to  allocate  s4\n"); 
exit  ( ) ; 
} 
} 
/*     */ 

/*   allocate  space  to  hold  columns  of  s4  for  passing  to  the  FFT  routine  */ 
/*     */ 

tempf ft= (COMPLEX*) calloc (p, sizeof (COMPLEX) )  ; 
if  (tempfft==NULL) { 

printf ("fam. ...  insufficient  space  to  allocate  tempfft   li  by  l\n",p); 
exit  ()  ; 
} 
/*    set  values  in  arguments  sent  to  fft   */ 
direction=l; 
norm=l; 

if  (cross==0)   {   /*  auto-fam   */ 
/*     */ 

/*   multiply  columns  of  s3   */ 
/*     */ 
for  (a=nprime-l;  a>  -1;  a — )  { 

for  (b=0;  b<nprime;  b++)   { 
/*   multiply  the  columns  of  s3  with  other  columns  of  s3    */ 
for  (j=0;  :<p;  j++) { 

tempfft [j] .r-(s3[j] [a] .r*s3[j] [b] ,r)+(s3[j] [a] .i*s3[j] [b] .i); 
tempfft [j] .i=(s3[j] [a] .i*s3[j] [b] .r)-(s3[j] [a] .r*s3[j] [b] .i); 
} 
/*  */ 

/*    FFT  EACH  COLUMN      */ 
/*     */ 

/*    go  get  FFT  performed  */ 
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f ft (tempf ft, p, direction, norm) ; 

/*  copy  tempf  ft  result  values  bacx:  into  a  row  of  s4  */ 

for  (j=(p/4);  j<((3*p/4)-l);  j++)  { 
c= ( (npnme-a-i) *p/2) +j- (p/4) ; 
s4[c] [b]=tempfft [j]; 
} 
}  /*  end  for  b=0. . .  */ 
}  /*  end  for  a=0 ...   */ 
}   /*  end  auto-fam   */ 
else  {   /*  cross-fam  */ 

*/ 
/*      multiply  columns  of  s3  with  y3   */ 
/*     */ 
for  (a=npnme-l;  a>  -1;  a — )  { 

for  (b=0;  b<nprime;  b++)   { 
/*   multiply  the  columns  of  s3  with  other  columns  of  y3    */ 
for  (]  =  0;  ]<p;  j++)  { 

tempfftt]]  .r=(s3[j]  [a]  .r*y3[j]  [b]  .r)  +  (s3[  j]  [a]  .i*y3  [  j]  [b]  .i) 
tempf  ft  [j]  .i-(s3[j]  [a]  .i*y3[j]  [b]  .r)-(s3[j]  [a]  .r*y3[j]  [b]  .i) 
} 
/*  V 

/*    FFT  EACH  COLUMN      */ 
/*     */ 

/*    go  get  FFT  performed  */ 

f ft (tempf ft, p, direction, norm) ; 

/*  copy  tempfft  result  values  back  into  a  row  of  s4  */ 

for  (j=(p/4);  j<((3*p/4)-l);  j++)  { 
c=( (nprime-a-1) *p/2) +j- (p/4) ; 
s4 [c] [b]=tempfft [j] ; 


}  /*  end  for  b=0 . . .  */ 
}  /*  end  for  a=0. . .   */ 
}   /*  end  cross-fam   */ 
/*     */ 
/*   OUTPUT   */ 
/*     */ 
halfp=p/2; 
if  (reduce==0)  {  /*  do  not  reduce  the  amount  of  output  */ 

if  ( ( (argc==8)&& (argv[7]  [l]=='a' ) )  I  I  ( (argc==7)&& (argv[6]  [l]=='a' ))) 
/*  make  an  ascii  output  file  of  the  whole  spectral  plane  */ 
/*  open  the  output  file  and  place  header  information  in  it  */ 
ofp  =  fopen (outf lie, "w" ) ; 
max_num_alf=npnme*p/2; 
max_num_f =npr lme ; 

fprintf  (ofp, "%i  %i\n",max_num_alf ,max_num_f )  ; 
for  (i=0;  i<max_num_alf ;  i++) { 
for  (j=0;  ]<nprime;  ]++) { 

fprintf (ofp, "%e   %e\n",s4 [i]  [j]  . r, s4[i]  [j]  . i)  ; 
}   /*  end  for  j=. . .   */ 
}   /*  end  for  i=. . .    */ 
f close (ofp) ; 
}   /*  end  if  ascii...   */ 

else  {  /*  make  a  plot  sxaf  no  reduced  output  file  */ 
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/*  open  i;he  output  file  and  place  header  information  in  it  */ 
ofp  =  fopen (outf ile, "w" ) ; 

/*  place  header  information  in  the  file  */ 
data_type=2;     /*  real_or_complex  */ 
max_num_alf=nprime*p/2;     /*  num_alpha  */ 
alpha_min=0 . 0;     /*  alpha_min  */ 
alpha_max=l . 0;     /*  alpha_max  */ 
max_num_f=nprime;     /*  max_num_f  */ 
f_min_all=  -0.5;     /*  f_mm_all  */ 
f_max_all=  0.5;     /*  f_max_all  */ 
fprintf  (ofp,  "%i  %i  %e  %e  %i  %e  %e\n",  data_type,  max_num_alf ,  alpha_mm, 

alpha_max,max_num_f , f_min_all, f_max_all)  ; 
for  (i=0;  Knprime;  i++)  {   /*  group  of  lines  */ 

for  (j=0;  ]<(p/2);  j++) {      /*    line  in  the  group   */ 
/*  place  alpha  line  header  information  in  file  */ 
curr_alf= (i*p/2) +];         /*  alpha  number  */ 
num_f=npnme-i;  /*  num_f   */ 

f_min=  -0  .  5+ (0  .  5*i/nprime)  ;  /*  f_mm   */ 
f_max=   0 . 5- (0 . 5*i/nprime) ;  /*  f_max   */ 

fprintf (ofp, "%i  %i  he   %e\n" , curr_alf , num_f , f_min, f_max) ; 
for  (k=i;  k<nprime;  k++) {    /*  points  on  the  line  */ 
tempi= (p/2* (nprime-1) )- ( (k-i) *p/2) +j; 
tempj=k; 

numreal=s4 [tempi] [temp]] .  r; 
numimag=s4 [tempi] [temp]] .i; 
fprintf (ofp, "he   %e\n", numreal, numimag) ; 
}   /*  end  for  k= . . .    */ 
}   /*  end  for  j=. . .    */ 
}   /*  end  for  i=. . .    */ 
/*   close  the  output  file   */ 
f close (ofp) ; 
}  /*  end  else  ...    */ 
}  /*  end  if  reduce=0  ...  */ 

else  {  /*  reduce  the  amount  of  output   */ 

if  ( ( (argc==9)&& (argv[7]  [i]=='a' ) )  I  I  ( (argc==8)&& (argv[6]  [l]=='a' ) ) )  { 
/*  make  an  ascii  output  file  of  the  whole  spectral  plane  */ 
/*  open  the  output  file  and  place  header  information  in  it  */ 
ofp  =  fopen (outf ile, "w") ; 
max_num_a 1 f =npr ime ; 
max_num_f=nprime; 

fprintf (ofp, "hi    %i\n", max_num_alf , max_num_f )  ; 
for  (i=0;  i<max_num_alf ;  i++)  { 
for  (j=0;  ]<nprime;  j++) { 
numreal=s4 [i*halfp] [j] .  r; 
numimag=s4 [i*halfp] [j] .  i; 

bigmag=sqrt (numreal*numreal  +  numimag* numimag) ; 
for  (k=l;  k<halfp;  k++)  {   /*  look  for  largest  value  */ 
tempreal=s4 [i*halfp+k] [j] .r; 
tempimag=s4 [i*halfp+k] [j] . i; 

tempmag=sqrt (tempreal*tempreal  +  tempimag*tempimag) ; 
if  (tempmag>bigmag)  {   /*  found  a  larger  value  */ 
b igmag= t empmag ; 
numreal=t empreal ; 
numimag=t  emp imag ; 
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}  /*  end  if  tempmag>bigmag. . .  *7 
}   /*  end  if  k=l  ...   */ 
fprintf (ofp, "he      %e\n", numreal, numimag)  ; 
}   /*  end  for  j=. . .   */ 
}   /*  end  for  i= . . .    */ 
f close (ofp) ; 
}   /*  end  if  ascn...   */ 

else  {  /*  make  a  plot_sxaf  reduced  output  file  */ 
/*  open  the  output  file  and  place  header  information  in  it  */ 
ofp  =  fopen (outfile, "w" ) ; 

/*  place  header  information  in  the  file  */ 
data_type=2;     /*  real_or_complex  */ 
max_num_alf=nprime;     /*  num_alpha  */ 
alpha_min=0 . 0;     /*  alpha_min  */ 
alpha_max=l . 0;     /*  alpha_max  */ 
max_num_f=nprime;     /*  max_num_f  */ 
f_min_all=  -0.5;     /*  f_mm_all  */ 
f_max_all=  0.5;     /*  f_max_all  */ 
fprintf  (ofp, "%i  %i  %e  %e  %i  %e  %e\n",  data_type,max_num_aif ,  alpha_mm, 

alpha_max, max_num_f , f_min_all, f_max_all)  ; 
for  (i=0;  i<nprime;  i++) {   /*  group  of  lines  */ 
/*  place  alpha  line  header  information  in  file  */ 
curr_alf=i;   /*  alpha  */ 
num_f=nprime-i;   /*  num_f   */ 
f_min=  -0  .  5+ (0  .  5*i/nprime)  ;    /*  f_mm   */ 
f_max=  0 . 5- (0 . 5*i/nprime) ;    /*  f_max   */ 
fprintf  (ofp, "%i  %i  he   %e\n", curr_alf , num_f , f_min, f_max) ; 
for  (k=i;  k<nprime;  k++) {    /*  points  on  the  line  */ 
tempi=  (halfp*  (npnme-1)  )  -  (  (k-i)  *halfp)  ; 
tempj=k; 

numreal=s4 [tempi] [tempj] .r; 
numimag=s4 [tempi] [tempj] . i; 

bigmag=sqrt (numreal *numreal  +  numimag * numimag ) ; 
for  (j=l;  j<halfp;  j++) {   /*  look  for  the  largest  value  */ 
tempreal=s4 [tempi+j ] [temp]] .  r; 
tempimag=s4 [tempi+j ] [temp]] . i; 

tempmag=sqrt (tempreal*tempreal  +  tempimag*tempimag) ; 
if  (tempmag>bigmag) {   /*  found  a  larger  value  */ 
b igmag=t empmag ; 
numreal=tempreal; 
numimag=t empimag ; 
}   /*  end  if  tempmag>  ....   */ 
}   /*  end  for  j=. . .    */ 
fprintf (ofp, "he   %e\n", numreal, numimag) ; 
}   /*  end  for  k=. . .    */ 
}   /*  end  for  i=. .  .    */ 
/*   close  the  output  file   */ 
f close (ofp) ; 
}  /*  end  else  binary...    */ 
}  /*  end  reduce  the  amount  of  output  */ 

}   /*  end  of  main  */ 
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APPENDIX  C.   SSCA  PROGRAM  USE 

Correct  commands  are: 

for  the  cross-ssca 

ssca  inputfilel  inputfile2  outputfile  n  nprime  1  Oflag 
ssca  inputfilel  input f ile2  outputfile  n  nprime  1  Oflag  red 

or  for  the  auto-ssca 

ssca  inputfilel  outputfile  n  nprime  1  Oflag 
ssca  inputfilel  outputfile  n  nprime  1  Oflag  red 

where:  inputfilel  is  the  file  containing  the  x  signal  samples 
inputfile2  is  the  file  containing  the  y  signal  samples 
outputfile  is  where  the  spectrum  values  will  be  placed 
n  is  the  number  of  samples  to  use  from  the  inputfiles 

and  it  must  be  a  power  of  2 
nprime  is  the  group  size  of  the  input  datasets  and  it 

must  be  a  power  of  2 
1  is  the  starting  offset  of  subsequent  datasets  and 

it  must  be  a  power  of  2 
Oflag  indicates  the  output  filetype,  ascii  or  binary 
-asc  for  an  ascii,  MATLAB  compatible  file, 

full  cyclic  spectral  plane 
-plo  for  a  plot_sxaf,    SSPI  compatible  file, 
half  cyclic  spectral  plane 
red  reduces  the  amount  of  data  output 

note  1:  the  first  two  entries  in  every  input  file  are 
expected  to  be  the  datatype  and  n.  datatype  =  1  for  real 
values,  2  for  complex  values.  n  is  the  number  of  samples 
contained  in  the  file,  one  sample  per  line. 

note  2:  error  may  occur  using  the  reduce  option  if 
(n-nprime) /nprime  is  not  an  integer. 

input  file  format: 

datatype    n 
sample  #1 
sample  #2 


sample  #n 
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output  ASCII  MATLAB  format 

nprime      n 
output (1)  (1) 
output  (1)  (2) 


output (nprime) (n) 


output  SSPI  plot_sxaf   format: 

datatype  (1  for  real,  2  for  complex) 
number  of  alphasalphaminalphamax 
number  of  f reqsf reqminf reqmax 

value  of  alpha  #1 

number  of  freqs  in  #lf reqmin_in_#lf reqmax_in_#l 

output  for  f reqmin_in_#l 


output  for  f reqmax_in_#l 


value  of  alpha  #alphamax 

number  of  freqs  in  #alphamaxf reqminf reqmax 

output  for  freqmin_in_#alphamax 


output  for  f reqmax_in_#alphamax 
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* include  <stdlib.h> 

^include  <stdio.h> 

#inciude  <math.h> 

^include  " /home3/carter/thesis/SSPI/pam/f ft . c" 

nnclude  " /home3/carter/thesis/SSPI/pam/ radix . c" 

main (argc, argv) 

int  argc; 

char  *argv[ ] ; 

{ 

COMPLEX  *x, **s3, **dwnconv; 

COMPLEX  *tempfft,  *sl,  *yl,  **s2,  *wmdow; 

int  i,  j,  k, type, n, nprime, p, 1, r, direction,  norm; 

mt  a,  cross,  file_n,m,  reduce,  sizes3,  redindex; 

int  tempi, temp], tempk, tempi; 

int  * out int; 

float  numreal,  numimag,  convf ac,  mamf ac; 

float  bigmag, tempmag, tempreal, tempimag; 

float  twopi=6. 28318530718; 

float  *outreal; 

double  z,y; 

char  *infilel,  *infile2,  *outfile; 

FILE  *ifpl, *ifp2, *ofp; 

/*   x  passes  data  to/from  fft  routine 

s3  holds  results  of  operations  after  first  fft 


si  holds  the 
yl  holds  the 


sample  values  from  mputfile  1 
sample  values  from  inputfile  2 


s2  holds  the  channelized  x  input  data  from  si 


check  for  the  correct  number  of  input  arguments 
the  correct  commands  are  : 

ssca  inputfile  outputfile  n  nprime  1  Oflag 
ssca  inputfile  outputfile  n  nprime  1  Oflag  reduce 
ssca  inputfilel  inputfile2  outputfile  n  nprime  1  Oflag 
ssca  inputfilel  inputfile2  outputfile  n  nprime  1  Oflag  reduce 
where:   inputfilel  is  where  the  x  signal  samples  are  expected  to  be 
mputfile2  is  where  the  y  signal  samples  are  expected  to  be 
outputfile  is  where  the  spectrum  values  will  be  put 
n  is  the  number  of  input  samples  to  use 
nprime  is  the  group  size  of  input  datasets 
1  is  the  starting  offset  of  subsequent  datasets 
Oflag  indicates  the  output  file  type 

-asc  for  an  ascn,  matlab  compatible  file 
-plo  for  an  ascn,  plot_sxaf  compatible  file 
reduce  is  an  option  to  generate  a  reduced  amount  of  output 

note  1:  type  and  file_n  are  expected  to  be  on  the  first  line  of  the  inputfile/s 
type  =  1  for  real,  2  for  complex, 
file_n  =  number  of  samples  following  the  first  line 

note  2:  error  may  occur  if  (n-nprime)  is  not  evenly  divisible  by  nprime 
when  using  the  reduce  option 


maintenance  history 

. . . .written  by  LCDR  Nancy 


Carter,  USN    NPGS  Monterey  CA   01  October  1992 
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.15  Oct  92 

.20  Oct  92 

.22  Oct  92 

.24  Nov  92 


.copied  from  current  version  of  fam.c  to  alter  to  ssca  algorithm 
modified  to  perform  ssca  per  Brown/Gardner/ Loomis  paper 
.made  output  file  ascn  format  MATLAB  compatible 
.made  an  output  format  SSPI  plot_sxaf  compatible 
.added  option  for  a  reduced  amount  of  output 


*/ 

/*  */ 

/*     CHECK  INPUT  VALUES     */ 
/*  */ 

if  (  (argc  !  =  7)&&(argc  !=  8)&&(argc  !=  9)  )  { 

print f ("ssca . . . fatal  error \n" ) ; 

printf (" incorrect  number  of  calling  arguments \n" ) ; 

printfC correct  formats  are:\n"); 

printf  (" ssca  mputfile  outputfile  n  nprime  1  Oflag\n"); 

printf (" ssca  input file  outputfile  n  nprime  1  Of lag  reduce \n"); 

printf (" ssca  mputfilel  inputfile2  n  outputfile  nprime  1  Oflag\n"); 

printf  (" ssca  mputfilel  mputfile2  n  outputfile  nprime  1  Oflag  reduce\n"); 

printf  (" wherexn"); 

printf  (" inputfilel  contains  the  x  signal  samples\n" ) ; 

printf  (" inputfile2  contains  the  y  signal  samples\n" ) ; 

printf  (" outputfile  will  contain  the  results\n"); 

printf  (" n  is  the  number  of  input  samples  to  use\n"); 

printf (" nprime  is  the  group  size  of  input  datasets\n" ) ; 

printf (" 1  is  the  offset  of  subsequent  datasets\n") ; 

printf (" Oflag  is  the  output  file  format\n"); 

printf  (" Oflag  =  -asc,  a  MATLAB  file  is  produced\n" ) ; 

printf  (" Oflag  =  -plo,  a  plot_sxaf  file  is  produced\n" ) ; 

printf  (" reduce  is  an  option  for  reduced  amount  of  output\n") ; 

exit  ()  ; 
} 
if  (argc==7)  {   /*  this  is  an  auto-ssca,  no  output  reduction   */ 

cross=0; 

mfilel=argv  [  1  ]  ; 

outfile=argv [2] ; 

n=atoi (argv[3] ) ; 

nprime=atoi (argv [ 4 ] ) ; 

l=atoi (argv [5] ) ; 

reduce=0; 
} 
if  ( (argc==8) && (argv(7] [0] !=' r' ) )  {   /*  this  is  a  cross-ssca,  no  data  reduction  */ 

cross=l; 

inf ilel=argv [ 1 ] ; 

mfile2=argv[2]  ; 

outfile=argv[3] ; 

n=atoi (argv [4] ) ; 

nprime=atoi (argv [5] ) ; 

l=atoi (argv [6] ) ; 

reduce=0; 
} 
if  ( (argc==8)&& (argv[7] [0]==' r' ) )  {   /*  this  is  an  auto-ssca,  with  data  reduction  */ 

cross=0; 

inf ilel=argv [1] ; 

outf ile=argv [2] ; 

n=atoi (argv [ 3 ] ) ; 

nprime=atoi (argv [4] ) ; 

l=atoi (argv [5] ) ; 
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reduce=l; 

} 

if  (argc==9)  {   /*  this  is  a  cross-ssca,  with  data  reduction  */ 

cross=l; 

mf  ilel=argv  [1  ]  ; 

inf ile2=argv [2 ] ; 

outf ile=argv[3] ; 

n=atoi (argv [ 4 ] ) ; 

nprime=atoi (argv [ 5 ] ) ; 

l=atoi (argv [ 6] ) ; 

reduce=l; 
} 

/*  verify  that  n,  nprime  and  1  are  powers  of  2 
i=n%2; 
if  (i!=0) { 

print f ("ssca . . . fatal  error\n" ) ; 

prmtf  (" calling  argument  n  is  not  a  power  of  2\n"); 

exit  ()  ; 
} 
i=nprime%2; 
if  (i!=0) { 

printf ("ssca . . . fatal  error\n" ) ; 

printf (" calling  argument  nprime  is  not  a  power  of  2\n") 

exit  ()  ; 
} 
j-1%2; 
if  (j!-0){ 

printf ("ssca. . . fatal  error \n") ; 

printf (" calling  argument  1  is  not  a  power  of  2\n"); 

exit  ()  ; 
} 
I*  */ 

/*   open  input  file  1,  prepare  to  get  the  x  signal  sample  data 
/*  */ 

lfpl  =  fopen (infilel, "r") ; 
fscanf (ifpl, "%i  %i", &type, &f ile_n) ; 
/*  verify  that  file_n  is  at  least  equal  to  n   */ 
if  (file_n  <  n) { 

printf  ("ssca . . . fatal  error \n" ) ; 

printf (" mputfilel  does  not  contain  enough  samples\n" ) ; 

fclose (ifpl) ; 

exit  ()  ; 
} 
/*  find  p  -  the  number  of  datasets (nprime  long)   */ 
p= ( (n-nprime) /l) ; 

/*  find  one  dimension  of  final  data  array  */ 
sizes3=p*l; 
/*     */ 

/*     READ  IN  DATA  SAMPLES    */ 
/*     V 

/*  allocate  space  to  read  in  the  x  sample  values  */ 
/*     */ 

sl= (COMPLEX*) calloc(n, sizeof (COMPLEX) ) ; 
if  (sl==NULL) { 

printf  ("ssca. .. insufficient  space  to  allocate  sl\n"); 

exit  ()  ; 
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if  (type==l)  {   /*  read  in  the  real  sample  values  */ 
for  (i=0;  i<n;  i++)  { 

fscanf  dfpl,  "%e\n",  inumreal)  ; 
si [i] . r=numreal; 
sl[i] .1=0.0; 
} 
} 

else  {   /*  read  in  the  complex  sample  values  */ 
for  (i=0;  i  <  n;  i++)  { 

fscanf (ifpl, "%e   %e\n", snumreal, &numimag) ; 
si [i] . r=numreal; 
si [i] . i=numimag; 
} 
} 

/*     */ 

/*   close  the  input  file  #1    */ 
/*     */ 
f close (ifpl) ; 
/*     */ 
if  (cross==l)  { 

/*  */ 

/*  open  input file2,  prepare  to  get  the  y  signal  sample  data   */ 
/*  */ 

ifp2  =  fopen (infile2, "r") ; 
fscanf (ifp2, "%i  %i", &type, &file_n) ; 
/*  verify  that  file_n  is  at  least  as  big  as  n    */ 
if  (file_n  <  n) { 

print f ("ssca. . . fatal  error\n") ; 

printf (" inputfile2  does  not  contain  enough  samples\n") 

fclose (ifp2) ; 
exit  ()  ; 
} 
/*     */ 

/*  allocate  space  to  read  in  the  y  sample  values  */ 
/*     */ 

yl= (COMPLEX* ) calloc (n, sizeof (COMPLEX) ) ; 
if  (yl==NULL) { 

printf ("ssca. .. insufficient  space  to  allocate  yl\n"); 
exit  ()  ; 
} 
/*  */ 

/*   read  in  the  complex  sample  values  */ 

/*  */ 

for  (i=0;  i  <  n;  i++)  { 

fscanf (ifp2, "%e   %e\n", snumreal, &numimag) ; 
yl [i] .r=numreal; 
yl [i] . i=numimag; 
} 
/*     */ 

/*   close  the  inputfile2    */ 
/*     */ 

fclose (ifp2) ; 
}   /*  end  if  cross  =1   */ 
/*     */ 
/*     FORM  DATASETS     */ 
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/*     */ 

/*  allocate  space  to  hold  the  p-by-nprime  ciatasets   */ 

*/ 
s2= (COMPLEX**) calloc (p,  sizeof (COMPLEX*) ) ; 
if  (s2==NULL) { 

printf ("ssca ...  insufficient  space  to  allocate  s2\n"); 
exit  () ; 
} 
for  (i=0;  i<p;  i++) { 

s2 [i]= (COMPLEX*) calloc (nprime, sizeof (COMPLEX) ) ; 
if  (s2[i]==NULL) { 

printf ("ssca ...  insufficient  space  to  allocate  s2\n"); 

exit  ()  ; 

} 

/*  */ 

/*  form  p  datasets  of  nprime  samples  from  the  original  sample  stream   */ 

/*  */ 

for  (i=0;  Kp;  i++)  { 

for  (j=0;  j<nprime;  j++) { 

k=i*l+j; 

s2[i] []]=sl[k]; 
} 
} 

/*  */ 

/*     APPLY  WINDOW  TO  DATASETS     */ 
/*     */ 

/*  allocate  space  to  hold  the  window  multiplicand  values  */ 
/*    */ 

window* (COMPLEX*) calloc (nprime, sizeof (COMPLEX) ) ; 
if  (window==NULL) { 

printf  ("ssca ...  insufficient  space  to  allocate  sl\n"); 
exit  ()  ; 
} 

/*  */ 

/*  calculate  the  window  values  for  the  "nprime"  sample  wide  window   */ 
/*  */ 

for  (i=0;  i<nprime;  i++) { 

y= (twopi*i) / (nprime-1) ; 

z=cos (y) ; 

window [i].r=  0.54  -  (0.46*z); 

window [i] .1=  window [i] .  r; 
} 

/*  */ 

/*  apply  the  window  to  rows  of  s2  containing  data       */ 
/*  */ 

for  (i=0;  Kp;  i++)  { 

for  (]=0;  ]<nprime;  j++) { 

s2  [i]  [  ] ]  . i=s2 [i]  [ j]  . i*window[ j]  . i; 
s2[i]  [j] .r=s2 [i]  [j] . r*window[ j] .r; 

} 
} 

/*  */ 

/*   FFT  EACH  DATASET  ROW     */ 
/*     */ 
/*   allocate  space  to  hold  rows  of  s2  for  passing  to  the  FFT  routine  */ 
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/*     */ 

x= (COMPLEX*) calloc (nprime, sizeof (COMPLEX) )  ; 
if  (x==NULL) { 

printf ( "ssca ...  insufficient  space  to  allocate  x\n"); 

exit  ()  ; 
} 

/*     */ 

/*   allocate  space  to  hold  rows  of  results  from  the  FFT  routine  */ 
/*     */ 

s3= (COMPLEX**) calloc (nprime, sizeof (COMPLEX*) )  ; 
if  (s3==NULL) { 

printf ("ssca. .. insufficient  space  to  allocate  s3\n"); 

exit  ()  ; 
} 
for  (i=0;  i<nprime;  i++) { 

s3[i]  =  (COMPLEX*) calloc ( (p*l) , sizeof (COMPLEX) )  ; 

if  (s3[i]==NULL) { 

printf ("ssca. .. insufficient  space  to  allocate  s3\n"); 
exit  ()  ; 
} 
} 

/*  */ 

/*  take  an  FFT  of  each  row  of  s2  */ 

/*  */ 

/*    set  values  in  arguments  sent  to  fft   */ 
direct ion=l; 
norm=l; 
for  (i=0;  Kp;  i++)  { 

/*    copy  row  of  s2  into  complex  array      */ 

for  ( j=0;  ]<nprime;  j++) { 
x[j]=s2[i]  [j]; 
> 

/*    go  get  FFT  performed  */ 

fft (x, nprime, direction, norm) ; 

/*  copy  x  values  into  1th  column  of  s3  */ 

r=i*l; 

for  (j=0;  j<nprime;  j++) { 
s3[j] [r]-x[j]; 

} 
} 

/*     */ 

/*     DOWNCONVERT  EACH  ROW   */ 
/*     */ 

/*   allocate  space  to  hold  the  downconversion  multiplicands   */ 
/*     */ 

dwnconv= (COMPLEX**) calloc (p, sizeof (COMPLEX*)  )  ; 
if  (dwnconv==NULL) { 

printf ("ssca. .. insufficient  space  to  allocate  dwnconv\n"); 

exit  ()  ; 
} 
for  (i=0;  i<p;  i++)  { 

dwnconv[i]= (COMPLEX*) calloc (nprime, sizeof (COMPLEX) ) ; 

if  (dwnconv[i]==NULL) { 

printf ("ssca. .. insufficient  space  to  allocate  dwnconv\n" ) ; 
exit  ()  ; 
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/*   downconvert  each  of  the  transform  sequences     */ 

*  *  / 

/*      calculate  the  down  conversion  multipliers     */ 
/  * 

mainfac=twopi*l/nprime; 
for  (m=0;  m<p;  m++)  { 
for  (k=0;  k<nprime;  k++) { 

convfac=mamfac*  (k-  (nprime/2)  +1)  *m; 
dwnconv [m] [k] . r=cos (convfac) ; 
dwnconv [m] [k] . i=  (-1 . 0) *sin (convfac) ; 
} 
} 

/*  multiply  downconversion  factor  against  each  frequency  value   */ 
for  (i=0;  i<p;  i++)  { 
tempi=i*l; 
for  (j=0;  ]<nprime;  ]++) { 

numreal=s3 [ j ]   [tempi]  . r* dwnconv [i]  [ j]  .  r  -  s3  [  j ]  [tempi]  . i*dwnconv [i]  [ j  ]  .  i; 
s3 [ j ]  [tempi]  . i=s3 [ j ]  [tempi]  . r*dwnconv [i]  [ j]  . i  +  s3  [  j ]  [tempi]  . i 'dwnconv [i]  [ j ]  . r; 
s3[j] [tempi] .r=numreal; 
} 
\ 

/*  */ 

/*  REPLICATE  COLUMNS     */ 

/*     */ 

for  (i=0;  i<p;  i++)  { 
r=i*l; 

/*  copy  column  into  1-1  adjacent  columns  */ 
for  (j=0;  ]<nprime;  j++)  { 
for  (k=l;  k<l;  k++)  { 

s3[j] [r+k]-s3[j] [r]; 
}   /*  end  for  k= . . .   */ 
}   /*  end  for  j=. . .   */ 
}   /*  end  for  i= . . .  */ 
/*     */ 

/*     MULTIPLY  COLUMNS     */ 
/*     */ 

if  (cross==0)  {   /*  auto-ssca   */ 
/*     */ 

/*   multiply  columns  of  s3  with  conjugate  value  of  si   */ 
/*     */ 

for  (a=0;  a<sizes3;  a++)  { 
for  (j=0;  ]<nprime;  ]++) { 

s3[j]  [a]  .r=(s3[j]  [a]  .r*sl[a]  .r)  +  (s3[  j]  [a]  .i*sl  [a]  .i) ; 
s3[j]  [a]  .i-(s3[j]  [a]  .i*sl[a]  ,r)-(s3[j]  [a]  .r*sl[a]  .i); 
} 
}  /*  end  for  a=0 ...  * / 
}      /*    end  auto-ssca  */ 
else  {   /*  cross-ssca   */ 
/*     */ 

/*   multiply  columns  of  s3  with  conjugate  value  of  yl   */ 
/*     */ 

for  (a=0;  a<sizes3;  a++)  { 
for  (j=0;  ]<nprime;  j++) { 

s3[j] [a] .r-(s3[j] [a] .r*yl[a] .r)+(s3[j] [a] .i*yl[a] .i); 
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s3[j] [a] .i=(s3[j] [a] .i*yl[a] .r)-(s3[j] [a] .r*yl[a] .i); 
} 
}  /*  end  for  a=0 ...   */ 
}   /*   end  cross-ssca   */ 

/*  */ 

/*    FFT  EACH  ROW     */ 
/*     */ 
/*     */ 

/*   allocate  space  to  hold  rows  of  s3  for  passing  to  the  FFT  routine  */ 
/*     */ 

tempfft= (COMPLEX* ) calloc (sizes3, sizeof (COMPLEX) ) ; 
if  (tempf ft==NULL) { 

printf ("ssca. .. insufficient  space  to  allocate  tempf ft   %i  by  l\n",sizes3) 
exit  ()  ; 
} 
/*    set  values  in  arguments  sent  to  fft   */ 
direction=l; 
norm=l; 
for  (j=0;  ]<nprime;  j++) { 

/*  copy  row  of  s3  into  tempfft   */ 
for  (k=0;  k<sizes3;  k++) { 
tempfft [k]=s3 [ j] [k] ; 
} 
/*    go  get  FFT  performed  */ 

fft (tempfft, sizes3, direction, norm) ; 

/*  copy  tempfft  result  values  back  into  a  row  of  s3  */ 

for  (k=0;  k<sizes3;  k++) { 
s3[ j] [k]=tempfft[k]; 
} 
}  /*  end  for  j=0. . .  */ 
/*  */ 

/*   OUTPUT  */ 

/*  V 

if  (reduce==0)  {  /*  no  output  reduction  */ 

if  ( ( (argc==7)&& (argv[6]  [l]=='a' ) )  I  I  (  (argc==6)&& (argv[5]  [l]=='a' ) ) )   { 
/*  make  an  ascn  output  file  full  spectral  plane  */ 
/*  open  the  output  file  and  place  header  information  in  it  */ 
ofp  =  fopen (outfile, "w") ; 
fprintf  (ofp, "%i  %i\n", nprime, sizes3)  ; 
/*   write  out  data  values  to  file  for  plotting   */ 
for  (i=0;  i<nprime;  i++) { 
for  (j  =  0;  :<sizes3;  j++) { 

fprintf  (ofp, "%e   %e\n",s3[i]  [j]  .r,s3[i]  [j]  .1) ; 
}   /*  for  j-0. . .  */ 
}   /*  for  i=0. . .  */ 
/*   close  the  output  file   */ 
f close (ofp) ; 
}  /*  if  ascn.  .  .  */ 
else  {   /*  make  a  plot_sxaf  output  file  half  spectral  plane  */ 
/*  open  the  output  file  and  place  header  information  in  it  */ 
ofp  =  fopen (outfile, "w" ) ; 
a=0;   /*  alpha  =  0  */ 

fprintf  (ofp, "%i  %i  %e  %e  %i  %e  %e\n", 2, sizes3, 0 . 0, 1 . 0, nprime, -0 . 5,  0 . 5)  ; 
for  (i=0;  i<nprime;  i++)  {   /*  group  of  lines  */ 

for  (j=0;  j< (sizes3/nprime) ;  j++)  {   /*  line  in  the  group   */ 
tempreal=  -0 . 5+ (0 . 5*i/nprime) ; 
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tempimag=  0 . 5- (0 . 5*i/npnme) ; 

fprintf (ofp, "%i  %i  %e  %e\n", a, (nprime-i) , tempreai, tempimag) ; 
a=a+l; 

for  (k=0;  k< (nprime-i) ;  k++) {    /*  points  en  the  line  */ 
tempi=npnme-k-l ; 
temp 3=  d+k)  *sizes3/nprime+]  ; 
numreal=s3 [tempi] [temp]] .r; 
numimag=s3 [tempi] [temp]] .  i; 
fprintf (ofp, "%e  %e\n", numreal, numimag) ; 
}  /*  end  for  k= . . .   */ 
}  /*  end  for  ]=. . .   */ 
}   /*  end  for  i= . . .   */ 
/*   close  the  output  file   */ 
f close (ofp) ; 
}  /*  end  if_else  */ 
}  /*  end  if  reduce=0...  */ 
if  (reduce==l)  {   /*  reduce  the  amount  of  output  data  */ 

if  ( ( (argc==8)&&(argv[6]  [ 1] ==' a' ) )  I  I  ( (argc==7)&& (argv[5]  [!]==' a' ) ))   { 
/*  make  an  ascn  output  file  full  spectral  plane  */ 
/*  open  the  output  file  and  place  header  information  m  it  */ 
ofp  =  fopen (outfile, "w") ; 
fprintf  (ofp, "%i  %i\n", nprime, p) ; 
redindex=l; 

/*   write  out  data  values  to  file  for  plotting   */ 
for  (i=0;  i<nprime;  i++) { 
for  (j-0;  ]<p;  j++) { 

numreal=s3 [i] [j*redindex] .r; 
numimag=s3  [i]  [j*redmdex]  . i; 

bigmag=sqrt (numreal *numreal  +  numimag * numimag ) ; 
for  (k=l;  k<redindex;  k++)  {   /*  look  for  largest  value  */ 
tempreal=s3 [i] []*redindex  +k] .r; 
tempimag=s3 [i] [ ]*redindex  +k].i; 

tempmag=sqrt (tempreal*tempreal  +  tempimag* tempimag)  ; 
if  (tempmag>bigmag)  {   /*  found  a  bigger  value,  swap  */ 
b igmag=t empmag ; 
numreal=tempreal; 
numimag=t emp imag ; 
}   /*  end  if  tempmag> . . . .  */ 
}  /*  end  for  k=l. . .   */ 
fprintf  (ofp, "%e   %e\n", numreal, numimag) ; 
}   /*  for  j-0. . .  */ 
}   /*  for  i=0. . .  */ 
/*   close  the  output  file   */ 
f close (ofp) ; 
}  /*  if  ascn ...  */ 
else  {   /*  make  a  plot_sxaf  output  file  with  reduction  half  spectral  plane  */ 
/*  open  the  output  file  and  place  header  information  in  it  */ 
ofp  =  fopen (outfile, "w") ; 
redindex=sizes 3 /nprime; 
a=0;   /*  alpha  =  0  */ 

fprintf (ofp, "%i  %i  %e  %e  %i  %e  %e\n", 2, nprime, 0 . 0, 1 . 0, nprime, -0  .  5,  0  .  5)  ; 
for  (i=0;  Knprime;  i++)  {   /*  nprime  alpha  levels  */ 
tempreal=  -0 . 5+ (0 . 5*i/nprime) ; 
tempimag=  0 . 5- (0 . 5*i/nprime) ; 

fprintf  (ofp, "%i   %i  he      %e\n", a,  (nprime-i) , tempreai, tempimag) ; 
a=a+l; 
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for  (k=0;  k< (nprime-i) ;  k++)  (    /*  points  at  the  alpha  level 
tempi=npnme-k-l ; 
temp  j=  d+k)  *redindex; 
numreal=s3 [tempi] [temp]] .  r; 
numimag=s3 [tempi] [temp]] .  i; 

bigmag=sqrt (numreal*numreal  +  numimag*numimag> ; 
for  (j  =  l;  ]<redindex;  j++)  {   /'*  line  in  the  group   */ 
tempreal=s3 [tempi] [temp]  +j] . r; 
tempimag=s3 [tempi] [temp]  +j] .  i; 

tempmag=sqrt (tempreal*tempreal  +  tempimag*tempimag) ; 
if  (tempmag>bigmag)  {   /*  found  a  bigger  value,  swap  */ 
b igmag=t empmag ; 
numreal=tempreal ; 
numimag=t  emp imag ; 
}   /*  end  if  tempmag>....  */ 
}  /*  end  for  ]=. . .   */ 

fprintf  (ofp, "%e  %e\n", numreal, numimag) ; 
}  /*  end  for  k=. . .   */ 
}   /*  end  for  i=. . .   */ 
/*   close  the  output  file   */ 
f close (ofp) ; 
}  /*  end  if_else  */ 
}  /*  end  if  reduce=l . . .  */ 
} 


fit 


APPENDIX  E.   SUBFAM  PROGRAM  USE 

Correct  command  is: 

subf am  n  inputfilel  offsetl  freqzerol  conj 1  inputfile2  offset2 
freqzero2  conj 2  outputfile  IOflag  Itypes 

where : 

n  is  the  number  of  samples  to  use  from  each  input  file, 
it  must  be  a  power  of  2 

inputfilel  is  the  first  file  containing  signal  samples 

offsetl  is  the  initial  sample  position  offset  location 
in  inputfilel 

freqzerol  is  the  flag  indicating  which  frequency 
modification  to  perform  on  inputfilel,  options  are: 
zn  -  zero  neg  freqs,  shift  pos  freqs  down 
zp  -  zero  pos  freqs,  shift  neg  freqs  up 
zz  -  don't  do  anything  to  frequency  components 

conj 1  is  the  flag  indicating  if  conjugation  of  the 
data  from  inputfilel  is  to  be  performed  prior  to 
multiplication  with  data  from  inputfile2,  options: 
y  -  yes,  conjugate  data 
n  -  do  not  conjugate  data 

inputfile2  is  the  second  file  containing  signal 
samples 

offset2  is  the  initial  sample  position  offset  location 
in  inputfile2 

freqzero2  is  the  flag  indicating  which  frequency 
modification  to  perform  on  inputfile2,  options  are: 
zn  -  zero  neg  freqs,  shift  pos  freqs  down 
zp  -  zero  pos  freqs,  shift  neg  freqs  up 
zz  -  don't  do  anything  to  frequency  components 

conj 2  is  the  flag  indicating  if  conjugation  of  the 
data  from  inputfile2  is  to  be  performed  prior  to 
multiplication  with  data  from  inputfilel,  options: 
y  -  yes,  conjugate  data 
n  -  do  not  conjugate  data 

outputfile  is  where  the  spectrum  values  will  be  placed 
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IOflag  indicates  the  datatypes  of  inputfilel, 

inputfile2  and  the  outputfile.  three  characters  are 
given  to  indicate  the  datatype  of  each  file. 
valid  datatypes  are: 

a  -  file  is  of  type  ascii 

b  -  file  is  of  type  binary 

ex.  "abb"  means  inputfilel  is  of  type  ascii, 

inputfile2  is  of  type  binary  and  the  outputfile 

is  to  be  of  type  binary 

Itypes  indicates  the  types  of  data  inside  inputfilel 
and  inputfile2.   the  outputfile  is  always  complex 
floating  point,  real  and  imaginary  components  on  the 
same  line.   two  letters  must  be  given  to  indicate 
the  datatype  for  each  input  file.   all  types  are 
converted  internally  to  complex  floating  point. 
valid  types  are: 
r  -  single  real  floating  point  samples  per  line 
i  -  single  integer  sample  per  line 
c  -  two  real  floating  point  numbers  per  line, 
representing  the  real  and  imaginary  parts 
of  a  single  sample 
ex.  "ic"  means  inputfilel  has  integer  samples 
and  inputfile2  has  complex  floating  point  samples 


example 


subfam  8192  datal.dat  64  zn  n  data2.dat  2156  zp  y 
outputA.dat  bbb  rr 

input  file  format:  output  file  format: 

sample  #1  magnitude  #1   phase  #1 

sample  #2  magnitude  #2   phase  #2 


sample  #n  magnitude  #n  phase  #n 
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APPENDIX  F.   SUBFAM  PROGRAM  LISTING 
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# include  <stdlib.h> 

^include  <stdio.h> 

^include  <math.h> 

4 include  "/home3/carter/thesis/SSPI/pam/f ft .  c" 

^include  " /home3/carter/thesis/SSPI/pam/ radix. c" 

/*  values  used  in  up/down  conversion  */ 

^define  INTFREQ  49375000.0     /*  intermediate  signal  freq  in  Hz  */ 

■[define  SAMPTIME  0.0000051282     /*  orig  signal  sampling  period  in  sec  */ 

main (argc, argv) 

mt  argc; 

char  *argv[] ; 

{ 

COMPLEX  *sl, *s2, -S3; 

int  i,  n,  direction,  norm,  of  f  1,  of  f2,  quarter  n,  scanmteger,  *readmteger,  *z; 

float  tempreal, tempimag, outputmag, outputphase, numreai, intfreq,  t,  convfac; 

float  *readreal, *readimag, *y; 

char  *mfilel,  *infile2,  *outfile,  *conjl,  *con]2,  *freqzerol,  *freqzero2,  *ioflag,  *itypes; 

FILE  *ifpl, *ifp2, *ofp; 

/* 

subfam  is  a  program  which  uses  signal  samples  from  two  input  files  to  calculate 
spectral  values  in  a  single  diamond  of  the  frequency-alpha  plane 

the  command  is  : 

subfam  n  mputfilel  offsetl  freqzerol  conjl  inputfile2  offset2  freqzero2 
conj2  outputfile  IOflag  Itypes 

n  is  the  number  of  samples  to  use  from  each  inputfile,  must  be  pwr  of  2 

mputfilel  is  the  name  of  the  first  input  file  containing  signal  samples 
offsetl  is  the  initial  sample  position  offset  location  in  mputfilel 
freqzerol  is  the  flag  indicating  which  frequency  modification  to  perform 
on  mputfilel,  valid  options  are: 

zn  -  zero  negative  frequency  components,  shift  positive 

components  down 
zp  -  zero  positive  frequency  components,  shift  negative 

components  up 
zz  -  don't  do  anything  to  frequency  components 
conjl  is  the  flag  indicating  if  conjugation  of  the  data  from  mputfilel  is 
to  be  performed  prior  to  multiplication  with  data  from  input file2. 
valid  options  are: 
y  -  yes  conjugate  data 
n  -  do  not  conjugate  data 

inputfile2  is  the  name  of  the  second  input  file  containing  signal  samples 
offset2  is  the  initial  sample  position  offset  location  in  mputfile2 
freqzero2  is  the  flag  indicating  which  frequency  modification  to  perform 
on  mputfile2,  valid  options  are: 

zn  -  zero  negative  frequency  components,  shift  positive 

components  down 
zp  -  zero  positive  frequency  components,  shift  negative 

components  up 
zz  -  don't  do  anything  to  frequency  components 
conj2  is  the  flag  indicating  if  conjugation  of  the  data  from  mputfile2  is 
to  be  performed  prior  to  multiplication  with  data  from  mputfilel. 
valid  options  are: 
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y  -  yes  conjugate  data 
n  -  do  not  conjugate  data 

outputfxle  is  the  filename  where  the  spectrum  values  will  be  placed 

lOflag  indicates  the  datatypes  of  input filel,  input file2  and  output  file. 

three  letters  must  be  given  to  indicate  the  datatypes  for  all  files 

valid  datatypes  are: 

a  -  file  is  of  type  ascn 

b  -  file  is  of  type  binary 

ex.  "abb"  means  mputfilel  is  of  type  ascn, 
inputfile2  is  of  type  binary,  and 
the  outputfile  is  to  be  of  type  binary 

Itypes  indicates  the  number  types  of  data  inside  mputfilel  and 

inputfile2.   the  outputfile  is  always  complex  floating  point, 

real  and  imaginary  components  on  the  same  line.   two  letters 

must  be  given  to  indicate  the  number  type  for  each  input  file. 

all  number  types  are  converted  to  complex  floating  point. 

valid  number  types  are: 

r  -  single  real  floating  point  sample  per  line 

i  -  single  integer  sample  per  line 

c  -  two  floating  point  numbers  per  sample  on  each  line 

ex.  "ic"  means  input filel  has  integer  samples  and 

inputfile2  has  complex  floating  point  samples 

sample  useage: 

subfam  8192  datal.dat  64  zn  n  data2.dat  2156  zp  y  outputA.dat  bbb  rr 

maintenance  history 
. . .  created  01  October  1992,  LCDR  Nancy  J.  Carter,  USN  at  NPGS  Monterey  CA 
...  10.15.92  -  added  freqzero  options  -  njc 
...  10.16.92  -  added  error  check  on  conjl  and  conj2  -  njc 
...  10.25.92  -  changed  down/up  shifting  to  same  code  as  in  old  version 


/*  check  for  the  correct  number  of  input  arguments  */ 

if  (argc  !=  13) { 

prmtf  ("subfam.  ..  fatal  error  ...  incorrect  number  of  calling  argument s\n")  ; 

printf ("\n") ; 

print f (".... correct  format  is:   \n"); 

printf ("   subfam  n  inputfilel  offsetl  freqzerol  conjl  mputfile2  offset2  \n"); 

printf ("        freqzero2  conj2  outputfile  IOflag  Itypes \n"); 

printf  ("\n") ; 

printf ("\n") ; 

printf (" n  is  the  number  of  samples  to  use  from  each  inputfile,  must  be  pwr  of  2\n" 

printf (" inputfilel  is  a  file  containing  the  signal  samples\n"); 

printf (" offsetl  is  the  initial  sample  position  offset  in  input f ilel\n" ) ; 

printf (" freqzerol  indicates  type  of  freq  zeroing  and  shifting  to  do  on  inputfilel\i 

printf  ("  zn zero  negative  freqs  and  shift  down\n"); 

printf ("  zp....zero  positive  freqs  and  shift  up\n"); 

printf ("  zz....do  not  zero  anything,  do  not  shift\n"); 

printf (" conjl  indicates  if  conjugation  of  data  from  inputfilel  is\n"); 

printf ("  desired  before  multiplication\n")  ; 

printf ("  y....yes  conjugate  data\n"); 
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print f("  n no  do  not  conjugate  data\n"); 

prmtf  ("  usual  setup  is  conjl  =  n,  conj2  =  y\n"); 

prmtf  (" input  file2  is  a  file  containing  the  signal  samples \n"); 

printf (" offset2  is  the  initial  sample  position  offset  in  input f ile2\n" ) ; 

printf (" freqzero2  indicates  type  of  freq  zeroing  and  shifting  to  do  on  inputfilel 

printf ("  zn zero  negative  freqs  and  shift  down\n"); 

printf ("  zp....zero  positive  freqs  and  shift  up\n"); 

printf  ("  zz....do  not  zero  anything,  do  not  shift\n"); 

printf (" conj2  indicates  if  conjugation  of  data  from  mputfile2  is\n"); 

printf  ("  desired  before  multiplication^" )  ; 

printf  ("  v.... yes  conjugate  data\n"); 

printf  ("  n no  do  not  conjugate  data\n"); 

printf ("  usual  setup  is  conjl  =  n,  conj2  =  y\n"); 

printf (" output file  is  where  the  spectrum  values  will  be  placed\n"); 

printf  (" lOflag  indicates  datatypes  of  mputfilel,  mputf ile2  and  output f  iie\n"  )  ; 

printf  ("  aaa ascii,  ascii,  ascii \n")  ; 

printf  ("  aab ascii,  ascii,  binary \n")  ; 

printf  ("  aba ascii,  binary,  ascii \n")  ; 

printf  ("  abb ascii,  binary,  binary \n"  )  ; 

printf  ("  baa ascn,ascii,ascn\n")  ; 

printf  ("  bab as  cii,  as  en,  binary  \n")  ; 

printf  ("  bba binary,  binary,  ascii \n "  )  ; 

printf  ("  bbb binary,  binary,  binary \n"  )  ; 

printf (" Itypes  indicates  number  types  of  input  samples\n" ) ; 

printf  ("  r single  real  float\n"); 

printf  ("  l single  integer\n"  )  ; 

printf  ("  c complex,  two  floats\n"); 

printf ("\n") ; 

printf (".... sample  useage : \n" ) ; 

printf ("  subfam  8192  datal.dat  64  zn  n  data2.dat  2156  zp  y  outputA.dat  bbb  rr\n"); 

exit  ()  ; 

} 
n=atoi (argv [1] ) ; 
mf  ilel=argv  [2  ]  ; 
of f l=atoi (argv [3] ) ; 
freqzerol=argv [4]  ; 
conjl=argv [5] ; 
mf  ile2=argv  [6]  ; 
of f2=atoi (argv [7] ) ; 
f reqzero2=argv [ 8 ] ; 
conj2=argv  [  9] ; 
outf ile=argv [10] ; 
lof lag=argv  [11] ; 
itypes=argv  [12 ] ; 

/*    verify   that   n    is   a   power   of   2  */ 

i=n%2; 

if    (i!=0)  { 

printf ("subfam. . . fatal  error\n") ; 

printf (" calling  argument  n  is  not  a  power  of  2\n"); 

exit  ()  ; 

} 
/*  verify  that  conjl  &  conj2  are  y  or  n   */ 
if  (  (  (conjl  [0]  !='y'  )&&  (conjl  [0]  !='n'  )  )  I  I  (  (conj2  [0]  !='y'  )&&(conj2[0]  !='n'  )  )  )  { 

printf ( "subfam. . . fatal  error\n" ) ; 

printf (" conjl  and  conj2  must  be  set  to  y  or  n\n"); 
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exit  ()  ; 
} 
/*  verify  that  ltypes  are  r,  i  or  c     *7 

if  (  (  (ltypes  [0]  !='  r'  )  &&  (ltypes  [0]  '.='  i'  )  &&  dtypes  [0]  !='c' ) )  I  I 
(  (ltypes  [1]  !='r')&&  (itypesfl]  !='  i'  )  &&  dtypes  [1]  !='  c'  )  )  )  { 
print f  (  "subfam.  .  .  fatal  error'xn"  )  ; 

pnntf  (" ltypes  must  be  set  to  r,  i  or  c\n"); 

exit  ()  ; 
} 
/*  */ 

/*  allocate  space  for  the  1-D  data  arrays   */ 
/*  */ 

si- (COMPLEX*) calloc (n, sizeof (COMPLEX) ) ; 
if  (sl==NULL) { 

pnntf  ("subfam.  ..  insufficient  space  to  allocate  sl\n"); 
exit  ()  ; 
} 
s2= (COMPLEX*) calloc (n, sizeof (COMPLEX) ) ; 
if  (s2==NULL) { 

pnntf  ("subfam.  ..  insufficient  space  to  allocate  s2\n"); 
exit  ()  ; 
} 
/*  allocate  space  to  read  in  binary  data  */ 

y= (f loat *)malloc (l*sizeof (float)  )  ; 
z=  (int*)malloc (l*sizeof (int ) )  ; 

/*  */ 

/*  INPUT  */ 

/*  */ 

/*  open  input  files,  get  the  signal  samples    */ 
/*  */ 

switch  (ioflag[0])  { 

case  'a'  :     /*  inputfilel  is  an  ASCII  file  */ 
if  ((  lfpl  =  fopendnfilel,  "r") )  ==  NULL  )   { 

pnntf  ("subfam.  .  .unable  to  open  ascn  inputfilel\n")  ; 
exit  ()  ; 
} 
/*   move  to  the  desired  starting  offset  position  in  the  file   */ 
for  (i=0;i<offl;i++) { 
switch  (ltypes [0] )  { 

case  ' r' :   /*  single  real  float  per  line  */ 
if  (fscanf  (ifpl, "%e", &tempreal) !=1) { 

pnntf  ("subfam.  .  .EOF  encountered  while  attempting  to  reach  offsetl\n") 
fclose (ifpl) ; 
exit  ()  ; 
} 

break; 
case  ' i' :   /*  single  integer  per  line  */ 
if  (fscanf (ifpl, "%i", scaninteger) ! =1) { 

pnntf  ("subfam.  .  .EOF  encountered  while  attempting  to  reach  offsetl\n") 
fclose (ifpl) ; 
exit  (); 
} 

breaks- 
default:   /*  complex,  two  floats  per  line  */ 

if  (fscanf (ifpl, "%e  %e", Stempreal, &tempimag) ! =2) { 

pnntf  ("subfam.  .  .EOF  encountered  while  attempting  to  reach  offsetl\n") 
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fclose(ifpl) ; 
exit  ()  ; 
} 
}   /*  end  switch  itypes[0]  */ 
}  /*  end  if  i=0;i<offl   */ 

/*  read  in  n  data  samples  */ 
for  (i=0; i<n; i++) { 
switch  (itypes[0])  { 

case  ' r' :   /*  single  real  float  per  line  */ 
if  (fscanf  dfpl,  "%e",&sl[i]  .r)  !=1)  { 

print f ("subf am. . .EOF  reached  while  reading  input f ilelXn" ) ; 
fclose (ifpl) ; 
exit  ()  ; 
} 
break; 
case  ' i' :   /*  single  integer  per  line  */ 
if  (fscanf  (ifpl, "%i", scan integer)  !=1) { 

printf ("subfam. . .EOF  reached  while  reading  input f ilelXn" ) ; 
fclose (ifpl) ; 
exit  ()  ; 
} 

sl[i] . r=scaninteger; 
break; 
default:   /*  complex,  two  floats  per  line  */ 

if  (fscanf  (ifpl, "%e  he", &sl [i]  .r, &sl [i]  . i) ! =2) { 

printf ("subfam. . .EOF  reached  while  reading  inputf ilelXn" ) ; 
fclose (ifpl) ; 
exit  ()  ; 
} 
}   /*  end  switch  itypes[0]  */ 
}   /*  end  for  i=0  i<n   */ 
fclose  (ifpl) ; 
break; 

case  'b':     /*  inputfilel  is  a  BINARY  file  */ 
if  ((  ifpl  =  fopendnfilel,  "rb")  )  ==  NULL  )   { 

printf ("subfam. . .unable  to  open  binary  input f ilel\n" ) ; 
exit  ()  ; 
} 
/*   move  to  the  desired  starting  offset  position  m  the  file   */ 
for  (i=0;  Koffl;  i++)  { 
switch  dtypes[0])  { 

case  ' r' :   /*  single  real  float  per  line  *7 
if  (fread(y,sizeof (*y) , 1, ifpl) 1=1) { 

printf ("subfam. . .EOF  encountered  while  attempting  to  reach  offsetlXn") 
fclose (ifpl) ; 
exit  ()  ; 

break; 
case  ' i' :   /*  single  integer  per  line  */ 
if  ( f read (z, sizeof (*z) , 1, ifpl) !=1) { 

printf ("subfam. . .EOF  encountered  while  attempting  to  reach  offsetlXn") 
fclose (ifpl) ; 
exit  ()  ; 
} 


TM 


Dec      9  12:23  1992   testalg/subfam. c  Page  6 


break; 
default:   /*  complex,  two  floats  per  line  */ 
if  (fread(y, sizeof (*y) , 1, if pi) 1=1) { 

print f ("subfam. . .EOF  encountered  while  attempting  to  reach  offsetl\n") 
f close  dfpl)  ; 
exit  ()  ; 


if  (fread(y, sizeof(*y),l,ifpl) !=1) { 

printf ("subfam. . .EOF  encountered  while  attempting  to  reach  offsetl\n") 
fclose  dfpl)  ; 
exit  ()  ; 
} 
/*   end  switch  itypes[0]  */ 
/*  end  if  i=0;i<offl     */ 


/*  read  in  n  data  samples  */ 
for  (i=0; i<n; i++) { 
switch  (ltypes  [0] )  { 

case  ' r'  :   /*  single  real  float  per  line  */ 
if  (fread(y, sizeof (*y) , 1, lfpl) !=1) { 

printf ("subfam. . .EOF  reached  while  reading  input f ilel\n" ) ; 
fclose  dfpl)  ; 
exit  ()  ; 
} 
si [i] .r=*y; 
break; 
case  ' i' :   /*  single  integer  per  line  */ 
if  (fread(z, sizeof (*z), 1, if pi) !=1) { 

printf ("subfam. . .EOF  reached  while  reading  input f ilel\n" ) ; 
fclose  dfpl)  ; 
exit  ()  ; 
} 

s  1  [ i ]  . r=  *  z ; 
break; 

default:   /*  complex,  two  floats  per  line  */ 
if  (fread(y, sizeof (*y), l,ifpl) !=1) { 

printf ("subfam. . .EOF  reached  while  reading  mputf ilel\n" ) ; 
fclose (ifpl) ; 
exit  ()  ; 
} 
sl[i] .r=*y; 
if  (fread(y, sizeof (*y), 1, ifpl) 1=1) { 

printf ("subfam. . .EOF  reached  while  reading  inputfilel\n") ; 
fclose (ifpl) ; 
exit  ()  ; 
} 
si [i] .r=*y; 
}   /*  end  switch  itypes[0]  */ 
}   /*  end  for  i=0  i<n   */ 
fclose (ifpl) ; 
break; 
default:      /*  invalid  format  specified  for  inputfilel  */ 

printf ("subfam. .. invalid  format  specified  for  mputf ilel\n") ; 
exit  ()  ; 
}  /*  end  switch  ioflag[0]  */ 
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switch  doflagtl])  { 

case  'a' :     /*  inputfile2  is  an  ASCII  file  */ 
if  ((  ifp2  =  fopen (infile2, "r") )  ==  NULL  )   { 

printf ("subfam. . .unable  to  open  ascii  input f ile2\n" ) ; 
exit  ()  ; 
} 
/*   move  to  the  desired  starting  offset  position  in  the  file   */ 
for  (i=0; i<off2; i++) { 
switch  (itypesfl])  { 

case  ' r' :   /*  single  real  float  per  line  */ 
if  (fscanf (ifp2, "%e", Stempreal)  !=1)  { 

printf ("subfam. . .EOF  encountered  while  attempting  to  reach  offset2\n") 
f close (ifp2) ; 
exit  (); 
} 

break; 
case  ' i' :   /*  single  integer  per  line  */ 
if  (fscanf (ifp2, "%i", scaninteger)  ! =1)  { 

printf ("subfam. . .EOF  encountered  while  attempting  to  reach  offset2\n") 
fclose (ifp2) ; 
exit  ()  ; 
} 

break; 
default:   /*  complex,  two  floats  per  line  */ 

if  (fscanf (ifp2, "%e  %e", Stempreal, &tempimag) !=2) { 

printf ("subfam. . .EOF  encountered  while  attempting  to  reach  offset2\n") 
fclose (ifp2) ; 
exit  () ; 
} 
}   /*  end  switch  itypes[l]  */ 
}   /*  end  if  i=0;Koff2   */ 

/*  read  in  n  data  samples  */ 
for  (i=0; i<n; i++) { 
switch  (itypesfl])  { 

case  ' r' :   /*  single  real  float  per  line  */ 
if  (fscanf (ifp2, "%e", &s2 [i] . r) ! =1) { 

printf ("subfam. . .EOF  reached  while  reading  input f ile2\n" ) ; 
fclose (ifp2) ; 
exit  ()  ; 
} 
break; 
case  '  i'  :   /*  single  integer  per  line  */ 
if  (fscanf (ifp2, "%i", scaninteger) !=1)  { 

printf ("subfam. . .EOF  reached  while  reading  inputf ile2\n" ) ; 
fclose (ifp2) ; 
exit  ()  ; 
} 

s2 [i] . r=scaninteger; 
break; 
default:   /*  complex,  two  floats  per  line  */ 

if  (fscanf (ifp2, "%e  %e", &s2  [i]  . r, &s2  [i]  . i)  ! =2) { 

printf ("subfam. . .EOF  reached  while  reading  inputf ile2\n" ) ; 
fclose (ifp2) ; 
exit  ()  ; 
} 
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}   /*  end  switch  itypes[l]  */ 

}   /*  end  for  i=0  i<n   */ 
f close (ifp2) ; 
break; 

case  'b' :      /*  input file2  is  a  BINARY  file  */ 
if  ((  ifp2  =  fopen(mfile2,  "rb")  )  ==  NULL  )   { 

printf ("subfam. . .unable  to  open  binary  input f ile2\n" ) ; 
exit  ()  ; 

/*   move  to  the  desired  starting  offset  position  in  the  file   */ 
for  (i=0;i<off2;i++) { 
switch  dtypes  [1]  )  { 

case  ' r' :   /*  single  real  float  per  line  */ 

if  (fread(readreal, sizeof (*readreal) ,  1,  ifp2)  !=1)  { 

printf ("subfam. . .EOF  encountered  while  attempting  to  reach  of f set2\n" ) ; 
f close (ifp2) ; 
exit  ()  ; 

break; 
case  ' i' :   /*  single  integer  per  line  */ 

if  (fread(readinteger,  sizeof  (*readmteger) ,  1,  ifp2)  !=1)  { 

printf ("subfam. . .EOF  encountered  while  attempting  to  reach  of fset2\n") ; 

f close (ifp2) ; 

exit  ()  ; 

} 

break; 
default:   /*  complex,  two  floats  per  line  */ 

if  (fread(readreal, sizeof (*readreal) , 1, ifp2) !=1) { 

printf ("subfam. . .EOF  encountered  while  attempting  to  reach  of fset2\n" ) ; 

fclose(ifp2) ; 

exit  (); 

} 
if  (fread(readimag, sizeof (*readimag) ,  1,  ifp2) ! =1)  { 

printf ("subfam. . .EOF  encountered  while  attempting  to  reach  offset2\n") ; 

f close (ifp2) ; 

exit  (); 

} 
}   /*  end  switch  itypes[l]  */ 
}   /*  end  if  i=0;i<off2    */ 

/*  read  in  n  data  samples  */ 
for  (i=0;i<n;i++) { 
switch  (itypes[l])  { 

case  '  r'  :   /*  single  real  float  per  line  */ 

if  (fread(readreal, sizeof (*readreal) , 1, ifp2) ! =1) ( 

printf ("subfam. . .EOF  reached  while  reading  mputfile2\n" ) ; 
fclose(ifp2) ; 
exit  ()  ; 
} 
s2 [i] . r=*readreal; 
break; 
case  '  i'  :   /*  single  integer  per  line  */ 

if  (fread(readmteger,  sizeof  (*readinteger) ,  1,  ifp2)  !  =1)  { 

printf ("subfam. . .EOF  reached  while  reading  input f ile2\n" ) ; 
f close (ifp2) ; 
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exit  () 


scaninteger=*readinteger; 

s2  [i]  .  r=scanmteger; 
break; 
default:   /*  complex,  two  floats  per  line  */ 

if  (fread(readreal, sizeof (*readreal) , 1,  ifp2)  !=1) { 

print f ("subf am. . .EOF  reached  while  reading  input f ile2\n" ) 
fclose (ifp2) ; 
exit  ()  ; 
} 
s2 [i] .r=*readreal; 
if  (fread(readimag, sizeof (*readimag)  ,  1,  ifp2)  !=1)  { 

printf ("subfam. . .EOF  reached  while  reading  mputf ile2\n" ) 
fclose (ifp2) ; 
exit  ()  ; 
} 
s2 [i] . i=*readimag; 
}   /*  end  switch  itypes[l]  */ 
}   /*  end  for  i=0  i<n   */ 
fclose (ifp2) ; 
break; 
default:      /*  invalid  format  specified  for  input f ile2  */ 
printf ("subfam. . .invalid  format  specified  for  input f ile2\n" ) ; 
exit  ()  ; 
}  /*  end  switch  ioflag[l]   */ 
/*  V 

/*  ZEROIZE  REQUESTED  FREQUENCY  SIDE  AND  SHIFT  APPROPRIATELY   */ 
/*  */ 

quarter_n=n/4; 
switch  (freqzerol [1] )  { 

case  ' n' :     /*  zero  negative  and  shift  down  */ 

direction=l;   /*  request  a  forward  transform  time-to-freq  */ 


norm=0; 

f ft (si,  n,  direction, norm) 

for  (i=0;  K(n/2)  ;  i++)  { 

sl[i] .r=0; 

sl[i] .1=0; 


/*  results  returned  in  si 
/*  zero  lower  half  */ 


direct ion=-l;   /*  request  an  inverse  transform  freq-to-time 
norm=2; 

f ft (si, n, direction, norm) ;  /*  results  returned  in  si  */ 
/*  downconvert  */ 

t=0.0;  /*  starting  time  */ 

for  (i=0;  i<n;  i++)  { 

convf ac=TWOPI * INTFREQ*t ; 

tempreal=cos (convfac) ; 

tempimag=  (-1)  *sm  (convfac)  ; 

numreal=sl [i]  . r*tempreal  -  si  [i]  . i*tempimag; 

si [i] . i=sl [i] . r*tempimag+sl [i] . i*tempreal; 

sl[i] .r=numreal; 

t=t+SAMPTIME; 


} 

break; 
case  'p' :     / 

direct ion=l; 


zero  positive  and  shift  up   */ 

/*  request  a  forward  transform  time-to-freq  */ 
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norm=0; 

f ft (si, n, direction, norm) ;   /*  results  returned  in  si   */ 

for  (i=(n/2);  i<n;  i++) {  /*  zero  upper  half  */ 

sl[i] .r=0; 

sl[i] .i=0; 

} 
direction=-l;   /*  request  an  inverse  transform  freq-to-time   */ 
norm=2; 

f ft (si,  n, direction, norm) ;  /*  results  returned  in  si  */ 
/*  upconvert  */ 

t=0.0;  /*  starting  time  */ 

for  (i=0;  i<n;  i++)  { 

convf ac=TWOP I  * INTFREQ*t ; 

tempreal=cos (convfac) ; 

tempimag=sm  (convfac)  ; 

numreal=sl [i] .r*tempreal  -  si [i] . i*tempimag; 

sl[i] . i=sl[i] . r*tempimag+sl [1] .i*tempreal; 

s 1 [ i ] . r=numreal ; 

t=t+SAMPTIME; 
} 

break; 
case  '  z' :     /*  do  nothing  */ 

break; 
default:      /*  invalid  format  specified  for  freqzerol  */ 

printf ("subfam. .. invalid  format  specified  for  f reqzerol\n" ) ; 
exit  ( ) ; 
}  /*  end  switch  */ 
switch  (freqzero2 [1] )  { 

case  'n' :     /*  zero  negative  and  shift  down  */ 

direct ion=l;   /*  request  a  forward  transform  time-to-freq  */ 

norm=0; 

f ft (s2, n, direction, norm) ;   /*  results  returned  in  s2   */ 

for  (i=0;  K(n/2);  i++)  {   /*  zero  lower  half  */ 

s2[i] .r=0; 

s2[i] .1=0; 

} 
direction=-l;   /*  request  an  inverse  transform  freq-to-time   */ 
norm=2; 

f ft (s2, n, direction, norm) ;  /*  results  returned  in  s2  */ 
/*  downconvert  */ 

t=0.0;  /*  starting  time  */ 

for  (i=0;  i<n;  i++)  { 

convf ac=TWOP I  * INTFREQ*t ; 

tempreal=cos (convfac) ; 

tempimag= (-1) *sin (convfac) ; 

numreal=s2 [i] . r*tempreal  -  s2  [i]  . i*tempimag; 

s2 [i] . i=s2 [i] .r*tempimag+s2 [i] . i*tempreal; 

s2 [i] .r=numreal; 

t=t+SAMPTIME; 
} 

break; 
case  'p' :     /*  zero  positive  and  shift  up   */ 

direction=l;   /*  request  a  forward  transform  time-to-freq  */ 

norm=0; 

f ft (s2,  n,  direction, norm) ;   /*  results  returned  in  s2   */ 

for  (i=(n/2);  i<n;  i++) {  /*  zero  upper  half  */ 
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s2 [l] .r=0; 
s2[i] .1=0; 
} 
direction=-l;   /*  request  an  inverse  transform  freq-to-time   */ 
norm=2; 

f ft (s2,  n, direction, norm) ;  /*  results  returned  in  s2  */ 
/*  upconvert  */ 

t=0.0;  /*  starting  time  */ 

for  (i=0;  i<n;  i++)  { 

convf ac=TWOP I  * INTFREQ*t ; 
tempreal=cos  (convfac) ; 
tempimag=sin (convfac) ; 

numreal=s2 [i] . r*tempreal  -  s2 [i] . i*tempimag; 
s2[i] . i=s2 [i] . r*tempimag+s2 [i] . i*tempreal; 
s2[i] .r=numreal; 
t=t+SAMPTIME; 
} 

break; 
case  '  z'  :     /*  do  nothing  */ 

break; 
default:      /*  invalid  format  specified  for  freqzero2  */ 

printf ("subfam. .. invalid  format  specified  for  freqzero2\n" ) ; 
exit  ()  ; 
}  /*  end  switch  */ 
/*        */ 
/*  MULTIPLICATION   */ 
/*  */ 

/*  create  space  for  the  output  array    */ 
/*        */ 

s3= (COMPLEX*) calloc(n,sizeof (COMPLEX) ) ; 
if  (s3==NULL) { 

printf ("subfam. .. insufficient  space  to  allocate  s3\n"); 
exit  ()  ; 
} 
if  ( (conjl [0]=='y' )&& (conj2 [0]=='y' ) )  {  /*  multiply  conjugate  si  times  conjugate  s2  */ 
for  (i=0;i<n;i++) { 

s3[i] .r=sl[i] .r*s2[i] .r-sl[i] .i*s2[i] .1; 
s3 [i] .i=(-l)*(sl[i] .i*s2[i] .i+sl[i] .i*s2[i] .r); 
} 
} 
if  ( (conjl [0] ==' y' )&& (conj2 [0]==' n' ) )  {  /*  multiply  conjugate  si  times  s2  */ 
for  (i=0; i<n; i++) { 

s3[i] .r=sl[i] .r*s2[i] .r+sl[i] .i*s2[i] .1; 
s3[i] .i=sl[i] .r*s2[i] .i-sl[i] .i*s2[i] .r; 
} 
} 
if  ( (conjl [0]==' n' )&& (conj2 [0]==' y' ) )  {  /*  multiply  si  times  conjugate  s2  */ 
for  (i=0; i<n; i++) { 

s3[i] .r=sl[i] .r*s2[i] .r+sl[i] .i*s2[i] .1; 
s3[i] .i=sl[i] .i*s2[i] .r-sl[i] .r*s2[i] .i; 
} 
} 
if  ( (conjl [0]=='n' )&&  (conj2 [0]=='n' ) )  {  /*  multiply  si  times  s2  */ 
for  (i=0;i<n;i++) { 

s3[i] .r=sl[i] .r*s2[i] .r-sl [1] . i*s2 [i] . i; 
s3[i] .i=sl[i] .r*s2[i] . i+sl [i] . i*s2 [i] . r; 


1DD 


Dec   9  12:23  1992   testalg/subfam. c  Page  12 


} 
/*  take  an  FFT  of  n-point  s3    */ 

*/ 
/*  set  values  in  arguments  sent  to  fft   */ 

direction=l;   /*  request  a  forward  transform  time-to-freq  */ 
norm=0; 

fft (s3,  n,  direction, norm) ;   /*  results  returned  in  s3  */ 
/*  */ 

I*    OUTPUT  */ 
/*  */ 

/*  take  magnitude  and  phase  of  s3  and  place  in  the  output  file  */ 
/*  */ 

switch  dof lag[2]  )  { 

case  '  a'  :     /*  outputfile  is  an  ascn  file  */ 
if  ((  ofp  =  fopen (out file, "w") )  ==  NULL  )   { 

print f  ("subfam.  .  .unable  to  open  ascn  outputf ile\n" ) ; 
exit  ()  ; 
} 
/*   write  out  data  to  the  ascn  file  */ 

for  (i=0;i<n;i++) { 

outputmag=sqrt ( ( s  3 [ i ]  . r *s3 [ i ]  . r )  +  (s3 [ i ] .i*s3[i]  . i ) ) ; 
outputphase=atan (s3[i].i/s3[i].r); 

fprintf (ofp, "%-16 . 6e  %-16. 6e\n", outputmag, outputphase)  ; 
} 
/*   close  the  output  file   */ 
f close  (ofp) ; 
break; 
case  'b':      /*  outputfile  is  a  binary  file  */ 
if  ((  ofp  =  fopen (outf lie, "wb") )  ==  NULL  )   { 

printf ("subfam. . .unable  to  open  binary  outputf ile\n" ) ; 
exit  ()  ; 
} 
/*   write  out  data  to  the  binary  file  */ 

for  (i=0;  i<n;  i++)  { 

outputmag=sqrt ( (s3 [i]  .r*s3 [i]  .r)  +  (s3[i] . i*s3  [i]  . i) ) ; 
outputphase=atan  (s3 [i]  . i/s3 [i]  .r); 
*y=outputmag; 
if  (f write (y,sizeof(*y),l,ofp) 1=1) { 

printf ("subfam. . .error  while  writing  outputf ile\n" ) ; 
f close (ofp) ; 
exit  ()  ; 
} 
*y=outputphase ; 
if  (f write (y,sizeof(*y),l,ofp) !=1) { 

printf ("subfam. .. error  while  writing  outputf ile\n" ) ; 
f close (ofp) ; 
exit  ()  ; 
} 
}   /*   end  if  i=0  ...   */ 
f close (ofp) ; 
break; 
default:      /*  invalid  format  specified  for  outputfile  */ 
printf ("subfam. .. invalid  format  specified  for  outputf ile\n" ) 
exit  ()  ; 
}   /*  end  switch  ioflag[2]   */ 
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*  end  of  mam  */ 


IDE 


LIST  OF  ABBREVIATIONS 

BPSK  binary  phase  shift  keying 

FAM  FFT  accumulation  method 

FFT  fast  fourier  transform 

FSM  frequency  smoothing  method 

PAM  pulse  amplitude  modulation 

QAM  quadrature  amplitude  modulation 

QPSK  quadrature  phase  shift  keying 

SCF  spectral  correlation  function 

SSCA  strip  spectral  correlation  algorithm 

SSPI  Statistical  Signal  Processing,  Inc. 

SUBFAM  sub- FFT  accumulation  method 
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